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Abstract 

A method  for  estimation  of  the  global  discretization  error  of 
solutions  of  operator  equations  is  presented.  Further  an  algorithm  for 
iterative  improvement  of  the  approximate  solution  of  such  problems  is 
given.  The  theoretical  foundation  for  the  algorithms  are  given  as  a 
number  of  theorems.  Several  classes  of  operator  equations  are  examined 
and  numerical  results  for  both  the  error  estimation  algorithm  and  the 
algorithm  for  iterative  improvement  are  given  for  some  classes  of  ordinary 
and  partial  differential  equations  and  integral  equations. 
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1 . Introduction 

A simple  technique  for  estimating  the  discretization  error,  or 
improving  the  accuracy  of  the  numerical  solution  of  a functional  equation 

F(y)  = 0 

is  discussed.  The  error  estimate  is  obtained  as  the  difference  between 
the  solution  of  the  discretized  problem 

<t>h(n)  = 0 

and  a slightly  perturbed  discrete  problem 

4>h  ( ) + $j^n)  = 0 ; 

the  improved  solutions  are  obtained  from  a sequence  of  perturbed  problems. 

The  perturbation  operator  is  a discrete  approximation  to  the 
operator  F that  is  of  higher  order  of  accuracy  than  (i.e.,  if  ^ is 
consistent  of  order  p with  F,  then  <j>^  is  consistent  of  order  q > p with  F, 
see  section  2 for  more  precise  statements).  We  can  view  -<j>^(n)  as  an 
approximation  to  the  local  discretization  error  of  the  operator  <t>^ . 

The  operator  ^ corresponds  to  a more  accurate  discretization 
method  than  4>^.  The  basic  requirement  on  <f>j^  is  the  accuracy;  the  operator 
does  not  even  need  to  be  stable  as  it  is  applied  only  to  the  known 
quantity  n to  produce  a constant  to  be  added  to  <t>^. 

The  perturbations  needed  to  iteratively  improve  the  solution  are 
obtained  from  operators  that  are  increasingly  accurate  approximations  to 
F,  and  the  negative  of  the  perturbation  can  again  be  viewed  as  increasingly 
accurate  approximations  to  the  local  discretization  error  of  the  operator 
fy,.  If  the  order  of  accuracy  of  the  operator  <j>h  is  p we  can  gain  p orders 
of  accuracy  per  iteration  if  the  perturbations  are  chosen  judiciously. 
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One  simple  and  direct  way  of  constructing  the  perturbations  is  the 
following  (other  ways  will  be  examined  in  section  4): 

For  the  error  estimation  algorithm  compute  the  residual  (on  a 
certain  set  of  discrete  points),  which  is  obtained  when  the  numerical 
solution  n of  the  discrete  problem  <j>^(n)  = 0 (a  solution  defined  only  on 
a set  of  discrete  points)  is  substituted  for  the  unknown  function  y in 
the  functional  equation.  Any  functionals  (e.g.  derivatives)  in  F(y)  = 0 
that  had  to  be  approximated  to  get  the  discrete  problem  4>^(n)  = 0 are 
approximated  by  linear  combinations  of  the  components  of  the  vector  n 
(solution  of  <f>h(n)  = 0)  when  the  residual  is  computed. 

To  improve  the  solution  the  same  technique  is  used,  but  now  the 
residuals  are  computed  as  the  sum  of  the  old  residuals  and  a new  residual 
computed  from  the  most  recent  approximation  to  the  solution.  Further 
increasingly  accurate  formulas  are  used  to  calculate  necessary  approxima- 
tions to  functionals  that  appear  in  F(y). 

This  technique  is  useful  if  there  exists  an  expansion  of  the 
global  discretization  error  (in  the  discretization  parameter  h)  of  the 
discretized  problem  <j>^(n)  = 0 to  sufficiently  high  order  and  with  smooth 
error  terms.  If  sufficient  smoothness  conditions  are  met  (e.g.  in  those 
cases  where  iterated  deferred  corrections  can  be  applied)  several  iterative 
improvements  can  be  made,  for  each  iteration  the  order  of  accuracy  of  the 
new  approximate  solution  is  increased,  in  the  same  way  as  when  iterated 
deferred  corrections  are  used. 

For  the  cases  where  one  cannot  perform  iterative  improvement  due  to 
the  existence  of  significant  non-smooth  terms  in  the  global  error  expansion 
one  may  still,  in  many  cases,  get  realistic  error  estimates  with  the 
technique  described. 
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Related  techniaues  are  difference  correction,  Fox  [1947],  Volkov 
[1957],  and  iterated  deferred  corrections  according  to  Pereyra.  Pereyra 
has  worked  extensively  on  boundary  value  problems  for  ordinary  differential 
equations,  Pereyra  [1967b,  1973]  and  has  produced  very  efficient  computer 
codes,  Lentini,  Pereyra  [1974,  1975a,  1975b]  for  such  problems.  He  has 
also  treated  linear  and  mildly  nonlinear  elliptic  boundary  value  problems, 
Pereyra  [1967b],  Pereyra  [1970]  and  analyzed  general  nonlinear  functional 
equations,  Pereyra  [1967a]. 

The  basic  computational  tools  for  high  order  approximations  to 
linear  functionals,  needed  both  in  deferred  corrections  and  iterative 
improvement  can  be  found  in  Ballester,  Pereyra  [1966],  Bjorck,  Pereyra 
[1970],  Galimberti,  Pereyra  [1970]. 

Whenever  iterated  deferred  corrections  can  be  used,  iterative 
improvement  can  also  easily  be  used.  However  to  be  able  to  perform 
iterated  deferred  corrections  one  must  know  and  be  able  to  approximate 
the  terms  of  the  local  error  expansion  for  the  discretization  operator 
<j>h,  while  for  iterative  improvement  that  is  not  necessary.  In  some 
cases,  e.g.  for  boundary  value  problems  y"  = f(x,  y) ; y(a)  = a,  y(b)  = 6, 
these  local  error  expansions  are  easily  found  (by  Taylor  expansions) 
and  involve  no  terms  that  are  difficult  to  approximate.  In  other  cases, 
e.g.  for  problems  where  F(y)  is  nonlinear  in  some  derivative  of  y,  the 
expansions  are  very  cumbersome  to  find  and  even  more  cumbersome  to 
approximate  (see  Pereyra  [1970]),  which  makes  iterated  deferred  corrections, 
although  theoretically  feasible,  in  practice  impossible  to  perform. 

Another  related  technique  is  iterated  defect  corrections  according 
to  Stetter  [1974],  Frank  [1975].  The  main  difference  between  the  approach 
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of  the  two  later  papers  and  our  approach  is  that  they  define  a smooth 
global  approximation  to  the  solution  of  the  operator  equation  F(y)  = 0 
from  the  discrete  solution  n of  ^(n)  = 0 and  then  compute  the  perturbation 
from  that  global  approximation,  while  we  use  the  discrete  solution 
directly  and  need  to  be  concerned  only  about  the  local  properties  of  the 
perturbations . 

In  section  3 of  Stetter  [1974]  the  author  also  introduces  a 
formalism,  similar  to  ours,  in  order  to  discuss  how  the  leading  term  of 
the  local  discretization  error  for  the  numerical  solution  of  ordinary 
differential  equations  usually  is  estimated.  He  does,  however,  not  use 
this  formalism  any  further. 

The  starting  point  for  our  investigation  was  Stetter  [1974]  and 
the  many  papers  of  Pereyra  (1967),  ...,  (1975)  and  our  primary  goal  was 
to  find  cheap  ways  to  estimate  the  errors  in  the  numerical  solution  of 
partial  differential  equations. 

Extensive  discussions  of  the  existence  of  asymptotic  expansions 
of  global  truncation  errors  for  the  numerical  solution  of  functional 
equations  can  be  found  in  Stetter  [1965],  Pereyra  [1967a],  and  Stetter 
[1973]. 


5 

2.  General  theory 
2.1  Preliminaries 

The  notation  of  this  section  is  greatly  influenced  by  Stetter 
[1965,  1973]  and  Pereyra  [1967a]. 

The  basis  for  the  results  discussed  here  is  the  existence  of 
asymptotic  error  expansions  for  the  global  discretization  error  of  the 
solution  to  the  discretized  problem.  For  a thorough  discussion  of  this 
matter,  see  the  references  above. 

We  consider  functional  equations 

(2.1)  F(y)  = 0 

where  F:  E -*■  is  a generally  nonlinear  operator  from  a Banach  space  E 
into  a Banach  space  E^1.  We  will  always  assume  that  (2.1)  has  a unique 
solution  y € E. 

For  the  purpose  of  numerical  solution  the  problem  (2.1)  is 
discretized  in  the  following  sense: 

We  define  families  - depending  on  a real  parameter  h 6 H, 

H = {h  = ^ , n integer  in  a subset  of  (v,  v 5 nQ}  with  nQ  > 0 fixed}  - of 

0 0 
Banach  spaces  E^,  E^J  and  of  bounded  linear  transformations  A^ , A^ , which 

map  E,  E°  into  Eh,  E^ , respectively. 

Then  we  choose  a family  of  (nonlinear)  operators  E^  -+  E^ 

such  that  for  z C E and  h 6 H or  for  z = y and  h € H 

(2.2)  *h(Ahz)  = Ah{F(z)  + 1 hV-yz)>  + 0(hM+1), 

v=p 

where  f : E -*■  E°  does  not  depend  upon  h. 

If  for  all  z £ E I | ( a^  z)  -a^F(z)||  Q = 0(hP)  the  operator  <t>h 

Eh 

is  said  to  be  consistent  of  order  p with  F. 
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The  expression  4>h(A^  y)  formed  with  the  solution  y of  (2.1)  is 
often  called  the  local  discretization  error  of 

The  original  problem  (2.1)  is  now  replaced  by  the  "algorithm" 

(2.3)  i>h(n)  = 0, 

which  is  supposed  to  have  a unique  solution  r\(h)  £ Eh  for  h € H.  See  the 
references  above  for  a discussion  of  the  unique  solvability  of  <t>h(n)  = 0. 
The  global  discretization  error  of  (2.3)  is  defined  as 

e(h)  = n(h)  - Ah  y 6 Eh 

where  y is  again  the  solution  of  (2.1). 

(2.3)  is  convergent  of  order  p (p  .>  1)  if 

l|e(h)|L  < C hP  for  h € H. 
n 

The  global  discretization  error  e(h)  admits  an  asymptotic  expansion 
to  the  order  M (M  ^ p)  if  there  are  ey  £ E,  v = p(l)M,  ev  independent  of  h, 
such  that 

M 

(2.5)  | | e ( h ) - Z h v * ey 1 1 ^ < CH  iv  for  h £ H. 

We  will  call  stable  at  A,  z if  there  exist  constants  S and  r > 0 such 
that  uniformly  for  all  h € H 

(2-6a)  lie'  - 52|Ie  s S[||*  (e')  - ♦ (e2)||  0 

h Eh 

for  all  $ , i * 1 , 2 such  that 

(2-6b)  ll*h(E1)  - $h(Ah  z) | | 0 < r 

Eh 

(cf.  def.  1.1.10  in  Stetter  [1973]). 
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To  estimate  the  global  truncation  error  choose  a family  of 
(non-linear)  operators  <f^ : Eh  -*•  such  that  for  z £ E and  h ( H 

(2.7)  $^(Ah  z)  = aJ{F(z)}  + 0(hq)  ; p < q * 2p 

When  the  solution  n of  problem  (2.3)  is  known,  solve  for  nE  from 

(2.8)  *h(nE)  + $E(n)  = 0 
Under  suitable  assumptions 

r q-1 

(2.9)  n - n = Ah[  z hv  ey]  + 0(hq) 

v=p 

In  this  report  we  are  mainly  interested  in  estimating  the  global 
truncation  error  for  given  algorithms,  however  the  idea  behind  (2.7), 

(2.8)  can  be  extended  to  iteratively  compute  better  approximations  to 
the  solution  of  (2.1 ). 

To  iteratively  compute  better  approximations  to  the  solution  of 
(2.1)  choose  a family  of  (non-linear)  operators  <t>^  ■ : E^  -*•  E®,  i =1,2, 

...  such  that  for  z ( E and  h £ H or  for  z = y and  h € H 

M 

(2.10)  .(a.  z)  = A°{F(z)  + Z hv  f (z)}  + 0(hM+1) 

M h h v=(i+l)p  7,1 

Put  n*^  = n and  compute  n1 , i = 1>  2,  ...  recursively  from 

(2.11)  ^(n1)  + Z * (nV_1 ) =0  i = 1,  2,  ... 

n v= 1 n,v 

Under  suitable  assumptions 

Mi  M.+l 

(2.121  ■Ahy*''h[v.(-„)phVe».f]*°(h'  ’ 
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A simple  example  of  the  use  of  the  operator  formalism  of  this 
section  is  given  in  note  5 after  theorem  2.  Readers  that  are  unfamiliar 
with  the  notation  of  this  section  are  advised  to  study  that  example  before 
they  read  the  theorems . 

N°te  1 ^(Ah  y)  = A°{F(y)}  + 0(hq) 

so  4>h  is  consistent  with  F of  order  q.  Further,  under  suitable  assumptions 
(see  the  theorems  below) 

*h(Ah  y)  + *h(n)  = Ah{F(y>}  + 0(hP) 

so  <p^  + 4>^(n)  is  consistent  with  F of  order  q.  We  can  view  as  an 

approximation,  accurate  to  order  q - 1 in  h,  to  the  local  discretization 
error  of  the  operator  <j>h> 

Note  2 *h,i(Ah  y)  = A°(F(y)}  + 0(h(i+1),p) 

so  . is  consistent  of  order  ( i + 1 ) • p with  F.  Further,  under  suitable 
assumptions,  (see  theorem  3 below) 

*h(Ah  y)  + \ = AS{p(y)}  + 0(h(i+1),p) 

(This  is  not  proved  in  the  theorem,  but  can  easily  be  proved  with  the 
technique  used  in  the  proof  of  the  theorem.)  Thus  <ph  + I <j>.  (nV”b 

n v=i  n’v 

is  consistent  of  order  (i+l)*p  with  F.  We  can  thus  view  - E <)>.  (nV~^) 

v=l  n,v 

as  an  approximation,  accurate  to  order  ( i +1 ) • p - 1 in  h,  to  the  local 
discretization  error  of  the  operator  <Ph. 


9 


Note  3 From  the  consistency  of  the  operators  above  and  the  stability  of 
4>h  (which  implies  the  stability  of  $h  + gh  for  any  constant  gh  £ eJ)  one 
can  derive  results  on  the  convergence  and  the  existence  of  asymptotic  error 
expansions  for  the  algorithms  (2.8)  and  (2.11). 

2.2  Basic  theorems 
Theorem  1 

Let  y be  the  unique  solution  of  F(y)  =v0  and  let 

a)  the  global  discretization  error  n - Ah  y of  (2.3)  have  an 

asymptotic  expansion 

Mo 

n - A.  y = a,{  E h^  e.}  + 6°(h) 
j=P  J 

o V1 

with  5 P + k and  ||<5u(h)||  = 0(b  u ) 

b)  the  expansions  (2.2)  and  (2.7)  hold  with  M >,  p and  p < q $ 2p 

c)  the  operator  <t>^  be  stable  at  Ah  y in  the  sense  of  (2.6) 

d)  there  exist  constants  L and  b such  that  uniformly  for  all  h £ H 

llfjly1)  - fj(y2)ll  o * Lily1  - y2llt 

for  all  y1  £ E,  i = 1 , 2 such  that 

I ly1  - yl lE  * b 

j = p,  p+l,  ....  N1  (N1  = min  (M,  MQ  - k,  q - 1) 

e)  there  exist  constants  C,  C^-  and  d such  that  uniformly  for  all 


h £ H 
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ll*hU')  - ♦„(«*)!!  o * c- Me1  - c2||e  -h-k 


- *hE(c2)||  0 * cE | |c ' - c‘||E  -h 

Eh  h 


1 ,2m  u-k 

El 


for  any  s1 , ?2  € Eh  such  that 


I U1  - \ y\ lE  5 d , i = 1,  2 


Then  the  solution  n of 


satisfies  the  inequality 


4>h  (nE ) + 4>E  (n ) = 0 


E 'V1 

nL  = A y + 0(h  1 ) 


where  = min  [M,  Mq  - k,  q - 1]. 
Proof 

Note  that 


0 = $h(nE)  + 4>^(n ) = 4>h ( nE ) - 4>h  (a^  y)  + 4>h ( Ah  y)  + 4>E ( n ) 

i .e. 

i’h (n^)  - 4>h(Ah  yX  = - <t>h(Ah  y)  ■ <^(0) 

C N.+l 

We  will  show  below  that  4>^(A^  y)  + = ) so  for  sufficiently 

small  h we  have 

I kh(nE)  - ^h(Ah  y) | | < r 
Thus  from  the  stability  of  <i>h  at  A^  y we  have 

MnE  - Ah  y 1 1 s S 1 1 4>h  ( nE ) - *h(Ah  y)||  = S||-  <t»h(Ah  y)  - <t>E(n)|| 

Here  and  in  the  sequel,  whenever  there  is  no  danger  of  confusion,  we  omit 
the  indices  on  the  norms  that  refer  to  the  actual  Banach  spaces. 
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Introduce  the  notation 

M°  j 

z = z hJ  e. 
j=P  J 

and  form 

♦ ^ (Ah  + ^h^  = + 

= 4>hUh  y)  - ^h(Ah(y  + z)  + 6°)  + *Jj(Ah(y  + z)  + 6°) 

f Mn+l-k 

= *h(Ah  y)  - *h(Ah(y + z))  + ^(Ah(y  + z))  + 0(h  ) 

Ni  . Ni 

= ^{^(y)  + z My)  hJ  - F(y  + z)  - r f.(y  + z) 

J=P  J j=P  J 

N.+l 

+ F(y  + z)>  + 0(h  1 ) 

o i Nl+1 

= A^  j:  CMy)  - My  + z)]  hJ>  + o(h  1 ) 

j=p  J 1 

N.+l 

= 0(h  1 ) 

Hence 

E V1 

n - Ah  y = 0(h  1 ) 

Corollary 

Let  the  conditions  of  theorem  1 be  satisfied  with  Mq  * 2p  + k - 1 , 
q = 2p,  M 5 2p  - 1,  then  the  solution  nE  of 

<t|h(nE)  + <t>^(  n)  = 0 

nE  = Ah  y + 0(h2p) . I 


satisfies  the  inequality 
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Theorem  2 


Let  y be  the  unique  solution  of  F(y)  = 0 and 

a)  conditions  a),  c)  and  e)  of  theorem  1 hold 

b)  F,  <(>h  and  ^ be  twice  continuously  Frechet  differentiable 

c)  the  inequalities  below  hold 


M 


>huh  y)  = A°{F(y)  + z f (y)  hj}  + 0(hM+1) 

j=p 


y)  = A°{F(y)}  + 0(hq)  P < q s 2p 

«-(Ah  y)  = A°{F'(y)}  + 0(hq’P) 

^‘(Ah  y)  = A°fF- (y) } + 0(hq-P) 

(2l  p (2) 

with  M 5 p.  Further  let  <i>^  '(a^  z)  and  4>h  (Ah  z)  be  uniformly  bounded 
with  respect  to  h for  any  z £ E.  Then  the  solution  nE  of 

<t>h(nE)  + = 0 


satisfies  the  inequality 

L N.+l 

n = Ah  y + 0(h  ) 

N-j  = min(M,  MQ  - k,  q - 1) 

Proof 

Use  the  same  notation  as  in  the  proof  of  theorem  1.  This  proof 
proceeds  in  the  same  way  as  that  proof,  except  that 
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<t>hUh  y)  + 4>^(n)  = $h(Ah  y)  - <J>h ( Ah (y  + z))  + <i>^(Ah (y  + z))  + 0(h 


V-'S 


*h(Ah  y)  - *h(Ah  y)  * *h(Ah  y)  Ah  z + *h(Ah  y) 

rl  N.+l 

+ 4>h  Uh  y)  Ah  z + 0(h  ) 


A°{-  F 1 (y)  z - ( E f;.(y)  hj)  z + F(y)  + F'(y)  z} 


J=P 


N.+l 

+ 0(h  1 ) 

N.+l 

= 0(h  1 ) 


Note  1 These  theorems  are  the  basis  for  estimation  of  global  discretization 

errors  for  algorithms  where  the  error  expansion  contains  only  a few  smooth 

terms.  If  k s 1,  and  p $ Mq<  p + k the  error  expansions  contain  at  least 

one  smooth  term  h^  e but  due  to  condition  e)  no  estimate  of  the  error  can 

P 

be  obtained.  Numerical  experiments  (see  section  4.5)  indicate  that  it  might 
be  possible  to  relax  this  condition.  No  theoretical  results  along  these 
lines  have  yet  been  obtained. 

Note  2 The  condition  on  M in  b)  of  theorem  1 is  somewhat  less  restrictive 
than  the  condition  on  Mg  in  a).  However,  to  obtain  the  asymptotic  expansion 
in  a)  one  would  need  to  have  M 5 p + k in  the  expansion  (2.2).  The  lower 
limit  on  M for  the  expansion  (2.7),  though,  may  be  of  value  in  some  cases. 
Note  3 The  constant  k of  condition  e)  is  in  general  the  order  of  the  highest 
derivative  present  in  the  functional  equation  F(y) . 

Note  4 These  theorems  differ  only  in  the  differentiability  assumptions  on 
the  operators  involved  and  the  set  of  elements  from  E for  which  the 
inequalities  (2.2)  and  (2.7)  must  be  valid. 
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Due  to  the  fact  that  (2.2)  and  (2.7)  only  are  needed  for  the  exact 
solution  y to  F(y)  = 0 in  theorem  2 some  very  interesting  types  of  perturba- 
tion operators  are  allowed  in  that  theorem,  while  they  are  not  allowed  in 
theorem  1.  In  essence,  all  operators  for  which  the  local  discretization 
error  is  of  order  greater  than  p can  be  used  in  theorem  2,  while  in  theorem 
1 the  order  of  consistency  of  4^  with  F must  exceed  p.  Several  examples  in 
section  4 will  clarify  this  distinction. 

Whether  the  more  relaxed  differentiability  conditions  of  theorem  1 
are  of  any  practical  importance  I don't  know.  Future  studies  will  hopefully 
provide  some  answers  to  this  question. 

Note  5 The  value  of  k of  condition  e)  in  theorem  1 depends  on  the  norm  for 
the  space  E®.  For  linear  multi  step  methods  for  initial  value  problems  for 
ordinary  differential  equations  one  can  (at  least  for  4^)  get  k = 0 by 
choosing  Spijker's  norm  (see  Stetter  (1973),  section  2.2.4,  p.  81-84)  for 
E®.  To  get  k = 0 also  for  4^  the  operator  must  be  chosen  judiciously,  cf. 
also  note  2 on  p.  41  of  section  4.1. 

Note  6 As  an  exercise  in  the  formalism  of  this  section  we  analyze  an 
algorithm  for  the  two-point  boundary  value  problem 

y"  = f(x,  y) 

y(a)  = a ; y'(b)  =6  ; f 6 C2s([a,  b]  x R)  ; b = b + e 

The  2s-continuity  of  f is  chosen  for  convenience  of  notation. 

1 . Operator  equation 

F:  E + E° 

with  E = C2s[a,  b] 

E°  = R * c[a,  b]  x R 


and  the  following  norms 
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M9llEo  = l9Vlg2l*a^|gWI 


g = g(x) , a s x s b € EL 


'z(a)  - a ’ 

F(z)  = z"(x)  - f(x,  z(x))  a s x s b 


z'(b)  - e 


2.  Discretization 


Introduce  the  grid  x.  = a + i h,  i = 0,  1 N + 1, 

h = 2(b  - a ) / ( 2N  + 1)  i.e. 


I I I 


I I I 

x0  X1  x2 


I I I L 


I I I I 


XN-1  XN 


and  discretize  the  problem  according  to 


«i+i  - 2*i + Vi 


- f(xi>  £.)  = 0 i = 1,  2 N 


- a = ° 


CN+1  " CN 


-6  = 0 


In  the  operator  formalism  this  means 


define 


V E-Eh  • Eh  = R 

Ah  Z = [z(xQ),  Z ( X i ) , ....  Z(XN+1)] 


E°  > eJ  ; E°  = R x RN  x R 
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for 


A 


Ah  g 


and  use  the  norms 


g(xi ) 


g(x) 

2 


1-1.2,  . 


a s x s b 


lull 


max 


OsuN+1 


for 


Then  define 


Ml  o ■ IV V lr2l  * ,^N 


y = 


i = 1,  2,  ....  N 


V Eh*Eh 


*h(E>  = 


£o  - » 


hn  - 2ti  * £i-i 


f(x1 , ci) ; 


N+l 
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Note  that,  by  Taylor  expansions  we  get 


z(x.+])  - 2z(x.)  + z(x.  .,) 


i-V 


- f(x.,  z(x. ) ) = z"(x.)  - f(x.,  z(x.)) 


. S;1  2 (2j+2)/  , , 2j 

(2j  + 2)!  z xi  h 


+ 0(h2s_1) 


and 


z(xMx1)  - z(xj  s-1 

— ^ - B - z'(b)  - b + ^ T2rHTT  h‘ 


wWl 


+ 0(h2s_1) 

for  z £ E.  If  we  define  f : E -*■  E®,  v = 2,  3,  2s  - 2, 


f (z)  = 0 
v' 


v odd 


V*>  * 


0 


( v V 2) ! z(v+2)(x)  a $ x ^ b 

2 _(v+l)/,.  wl\V+l 

(v'  + T)!  2 (b)(2} 


v even 


we  have 


(a,  z)  = a°{F(z)  + l f (z)  hv}  + 0(h2s_1) 


h 


v=2 


3.  Perturbation 


Define 
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i =1,2,  ....  N 


V Eh-R 


50?0  - 75^  - 20^2  + 70C3  - 30^  + 5C& 


for  1 = 1 


DiU)  » 


-S^.2  + 80S...J  - 150Ci  + 80?i+1  - 5^i+2 


for  i = 2,  . . . , N - 1 


5^N-4  ' 30CN-3  + 70^N-2  " 20CN-1  " 75?N  + 50CN+1 


for  i = N 


and  Rk:  E.  -*■  R 

b h 


Rb(0  ■ 


'5N-3  + 5^N-2  ‘ 95N-1  " 17CN  + 22CN+1 


Note  that,  by  Taylor  expansions  we  get 


D^z)  = z"(x.)  + 0(h4)  i = 1,  2, 


Rb(Ah  z)  = z 1 ( b)  + 0(h*) 


for  z £ E . 


Now  define 


- « 


;U)  = MO  - f(x,,  t.)  i = 1,  2 N 


then  for  z £ E 


MO  - B 


$E(a.  z)  = A°{F(z)J  + 0(h4) 


4.  As  a further  exercise  we  examine  some  of  the  conditions  of  theorem  1 


f I 
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"Vy  > - Vy  0 


„'(W)  . 

n +2iyr  (y1(jtl)(b)  - y2(Jt,,(b))-( 


+ 2 ™*_  ly1(j*g)(«)  - 

asxsb  ( j + 2) ! 

s 2-(]-)j+1  lly1  - y2 1 1 E + 2-|  ly1  - y2| 

* 3* | ly1  - y2|  L 


,2( j+2) , 


J .2 


for  any  y i y ( E,  j = 2,  4,  ....  2s  - 2. 


Further 


My1)  - My2) II  n = o 


Lj+i 

2 


for  j = 3,  5,  2s  - 3 

so  condition  d)  of  theorem  1 is  satisfied. 
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m 


i 


- ‘h^HI  0 


c>  - E2 


‘w - «f)  ♦«!.,- «i.i  ...  , 

T (f(x,,  £,. 


1 ’ ’1 


C1  - C2  - (C1  - £2) 
^N+l  ^N+l  UN 


_ iJ  j-2,  A !cN+1  ' 5N+1  " (CN  ‘ $ 
'*0  ' ^O1  h 


c1  - £2  - 2ft1  - r2l  + r1  - E2 

+ max  [ i + 1 S'+l  Ci-1  ^i-1 


1 <h  N 


- f (x. 

y i 


lie1  - e2||E  + £ lie1 


< c- 1 u1  - £2||  -h' 

n 

where 


C • h‘(l  * My)  * 2h0  * 4 


£,  € intU],  £•] 


M = m3X 

y a j x $ b 
|z  - y(x) | * d 

Analogously  we  get 


f (x,  z(x) | ; lie1  - A y | | : 

y n th 


ll*h(e1)  - *2(c2) 


s Cf-llc'  - £2Mc  •!>' 


where 


\ 

)-f(x(.t2)  II  0 

n 

/ 

. ^)(e]  - e2)| 
lie1  - e2| |E 

th 


/ I 
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C 


E 


*0<’ 


M ) 

y 


54 

24 


320 

60 


1 2 

and  M , q , ^ are  as  above.  Thus  condition  e)  is  satisfied. 

Theorem  3 

Let  y be  the  unique  solution  of  F(y)  =0  and 

a)  the  global  discretization  error  n ■ y of  (2.3)  have  an 
asymptotic  expansion 


n 


Ah  y = A.{  z hj  e.}  + 6°(h) 
j=P  J 


with  M s u(p  + k)  and  J | 6° | | = 0(hM+1) 

b)  the  expansions  (2.2)  and  (2.10)  hold  with  M 5 u(p  + k) 

c)  the  operator  <f>h  be  stable  at  Ah  y in  the  sense  of  (2.6) 

d)  the  operators  F,  f.  and  f.  . have  the  following  differentiability 

J > > J 

properties : 

F:  [M/p]  + 1 times  Frechet  differentiable 

f^:  [(M  - j)/p]  + 1 times  Frechet  differentiable 

J 

j=p,p+l,...,M 

fVj:  C(M  - j)/(vp)]  + 1 times  Frechet  differentiable 

j = (v  + 1 ) • p , (v  + 1 ) • p + 1,  ....  My 
v = 1 , 2,  . . . , p - 1 

Here  [expression]  stands  for  the  integer  part  of  the  expression  within  the 
square  bracket. 

e)  there  exist  constants  d,  C and  C , v = 1 , 2,  . . . , u - 1 such 
that  uniformly  for  all  h € H 
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I Uh (C1)  - *hU2)||  < ells1  - C2||-h'k 
' ^h,v(c2)|1  f Cvllc1  ‘ S2H*h'k 


for  all  £ Eh,  1 = 1 , 2 such  that 


I U1  - Ah  y|  |E  s d 


f)  [F'(y)]"1  exist 

g)  it  be  possible  to  define  e j = (v  + 1)  p M , v = 1,  2, 

v j V 

...,u-l,M=M-k-v  according  to 

F 1 (y)  e . = q . 
vj  avj 

where  g^  are  defined  in  equation  (2.18)  below;  then  the  global  discretiza- 
tion error  n1  - Ah  y of  (2.11)  has  an  asymptotic  expansion  of  the  form 

M. 

n’  • Ah  y = V . ,?  eii  hJ}  + 6’(h) 


j=(i+l ) *p  1J 


M.+l 


where  | ( 6 1 1 1 = 0( h 1 ) ; M.  = M - ik,  i = 1 , 2 y - 1 

Proof 

Introduce  the  family  of  operators  ^ . : Eh  -*•  E®,  i * 1 , 2,  , 

y - 1 by 

(2.13)  *h,1(t)  - *„<C)  * I *h,j(nj'') 

J ’ 

Then  n1 , i = 1 , 2,  . . . , y - 1 are  the  solutions  of 

(2.14)  i|>h  j (n^ ) = 0 
Note  that  as 
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0 = Vi-l^1"1)  = V*1"1)  + lj1  % 

j=l 

we  have 

0 = ^h,i(nl  > = *h(r,1)  ' *h(n1_1)  + *h,i(rii"1) 


0h(ri1)  = $h(n1"1)  - ^hji(n1_1) 

Further  note  that  as  and  only  differ  by  an  additive  constant  the 
stability  of  is  implied  by  the  stability  of  ify 

We  will  prove  the  theorem  by  induction  on  i,  so  assume  that 


(2.15) 


i-l  . 

Ah  y = Ah{.  Z ei-l  j hJ}  + 6 (h) 

j=i*p  J 


i 1 M.  +1 

with  | |6  ||  - 0(h  ) and  M.  = M - i-k. 

Introduce  the  notation 


(2.16) 


zs  = I e . h^ 


S=0,  1,  . . . , y - 1 


si  * '»  •••*  h ■ 

J= (S+1 ) • p SJ 

From  (2.14),  (2.13),  (2.2)  and  (2.10)  we  get 

l K,i(n1)  • Vi(Vy  + z1})H  = 0(hP) 

so  for  sufficiently  small  h we  have 

• ' ^h,i(Ah{y  + z1})N  < r 

Hence  from  the  stability  of  . we  get  for  h £ H 

n » i 

ll6l||  = 1 1 n 1 - Ah  z1  - Ah  y|  | = | In1  - iyy  + z 1 } | | 

* sl  kbjtn1)  - ,i  (Ah(y  + z1})  1 1 
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M.+l 


To  prove  that  Hs1!!  = 0(h  1 ) if  e.  . are  defined  as  in  g)  we  note  that 

• fj 

Vi(r,1)  ' *h,i(Vy  + z1})  = *h(n1)  - Wy  + z'}) 

= 4>h(n3  b - 4>h ^i (n1  -1 ) - ( Ah {y  + zi}^ 


M,  ,+1-k 


* *h<Vy  + z1"  })  ’ *h,i(Vy  + z1'1})  ' *h(Ah{y  + z1})  + 0(h  i_1  } 

M. 

= A°{F(y  + z1'1)  + i1  f (y  + zi_1)  hj  - F (y  + z1"1) 


J=P 


Mi  • , Mi  M.+l 

l f-  ,(y  + z ) h3  - F(y  + z1)  - z f . (y  + z1)  hJ}  + 0(h  1 ) 

j=(i+l)*P  1,J  j=p  3 

M. 

!(-  F(y  + z1)  + z [f.(y  + z1"1)  - f (y  + z1’)]  hj 
j=P  J J 


M. 


1 ill  V1 

I Ay  + Z1  ')  hJ)  + 0(h  1 ) 

j=p(i+l)  1,3 


n.  M. 

- A°{F(y)  + z1  JTF(r,(y)(zi)r  - z ( z ^ f<r) (y)[(zi_1  )r  - ( z1  )r])  hj 
r= 1 r-  j=p  r=l  r>  3 


M. 


+ i1  [ z ^ f jrhy)(zi"1)r]  h3}  + 0(h  1 ) 

j=p(i+l ) r=0  r’  1,3 

The  upper  summation  limits  in  the  two  sums  over  r above  are  chosen 
such  that  all  relevant  terms  are  included. 

Insert  the  expressions  for  z1  and  z1  and  collect  terms  of  equal 
powers  of  h,  then 

Mi 

*h,1(n'>  ‘ *M(V*  * z'’>  = - 1h<j  = (?*1)p(F'<y)  6i  j ' 9i  j>  hJ> 
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where  g.  . is  independent  of  e.  . and  h and  can  be  constructed  from 


J 


i J 


(2.18) 


M.  M.  M. 

z g,-,-  hj  = - D1  pr  F(rJ(y)(  / e.„  h*)r 


j=(i+l)p  1J  r=2 


2=(i+l )p 


'u 


M. 


Mi-1 


+ S ( S 7T^r)(y)[(  1 ' ei-U  h^r 
m=p  r=l  r>  m 2=ip  1 11 


M. 

l 


i e.  h£)r])hm 


t=(i+l)p 


12. 


- i [ 2 ]t  fjl)(y)(  i 1 e1.u  hV].hm 

m=(i+l )*p  r=0  r-  im  2=i p 1 n 

♦ o(hMi+1) 


Now,  if  e.  . , j = (i+l)p,  ....  M.  are  the  solutions  of 

• J 


- 9^-0 


then 


^h,i(n1)  - *h,1(Vy  + z1})  = 0(h  1+1) 

and  thus  from  (2.17) 


1 1*1 1 1 = o(hMi+1) 

But  from  condition  a)  the  assumption  of  induction  is  true  for  i = 0 with 
= n,  thus  the  theorem  is  proved.  Q.E.D. 

Theorem  4 

Let  y be  the  unique  solution  of  F(y)  = 0 and  let 

a)  all  the  conditions,  except  b),  of  theorem  3 hold 

b)  $h,  <t^  be  [M/p]  + 1 times  continuously  Frechet  differentiable 

c)  the  inequalities  below  hold  for  v = 0,  1,  ....  [M/p] 
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M 0 M M"v*p  ful  i rn-i-vp 

K y)  = }(y)  + £ nv,(y)  hJ}  + 0(h  ) 


M+l-vp 


'h  vuh 


J=P 

M-vp 


♦iv}(Ah  y)  = A^p(v)(y)  + 1 fiV!(y)  hJ> + o(hM+1-vp) 

h>1  h h j=(i+i).p  1,j 

for  i = 1 , 2,  . . . , p - 1 

then  the  global  discretization  error  n1  - y of  (2.11)  has  an  asymptotic 
expansion  of  the  form 

M. 

i 1ii 

n ■ a.  y = A.  { z e. . hJ}  + ^(h) 

h hj=(i+l)p  1J 


where 


1 1*1 1 1 * 0(h  i+1)  ; M.  = M - i*k 

i = 1,  2,  ....  w - 1 

The  proof  of  this  theorem  is  analogous  to  the  proof  of  theorem  3,  cf.  the 
proofs  of  theorem  1 and  2. 


Note  1 These  theorems  are  the  basis  for  iterative  improvement  of  the 
numerical  solution  of  an  operator  equation.  When  the  conditions  of  the 
theorem  are  satisfied  at  least  u - 1 improvements  can  be  made.  At  no 
extra  expense  an  estimate  of  the  global  discretization  error  of  the 
solution  n1  can  be  obtained  as  n^  - n1-^. 

Note  2 The  condition  on  M in  the  expansion  (2.10)  can  be  relaxed 

somewhat,  and  could  be  made  dependent  on  i so  M would  decrease  with  i. 

combinations  of  the  conditions  on  <t>,  and  <j>.  of  theorems 

h h,v 

d used.  We  can,  e.g.,  use  the  conditions  of  b)  and  c)  on 

i>om  theorem  4 and  the  conditions  on  <t>,  of  d)  from  theorem  3. 
n , v n 
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Note  4 The  main  problem  of  applying  this  theorem  is  of  course  the 
construction  of  the  operators  r.  Several  examples  of  such  operators 
will  be  given  in  section  4.  To  illustrate  a major  difficulty  and  resolve 
it  we  will  again  consider  the  two  point  boundary  value  problem 

y"  = f(x,  y) 


y(a)  = a y 1 (b)  = s 

of  note  5 to  theorem  2.  We  cannot  use  the  operator  as 

D1(Ah  z)  = z" ( x i ) + c^z^^x^-h4 
Di(Ah  z)  = z"(x.)  + c0*z^6^(x.)*h4 
DN(Ah  z)  = z"(xN)  + c2-z(6)(xN)-h4 


^ because 


C1  * c0  * c2 


so  the  expansion  (2.10)  is  not  valid.  This  is  caused  by  the  use  of 
approximation  formulas  for  the  two  points  closest  to  the  boundary 
that  are  different  from  the  formulas  in  the  interior  of  the  interval. 
However,  if  we  fix  the  number  of  iterative  improvements  that  we  want  to 
make  to  (u  - 1)  we  can  construct  operators  D.:  + R,  i s 1,  2 N 


DjU) 


N+l 

I 

v=0 


TV 


V 


(some  a.  v may  be  zero)  such  that 

D.(Ah  z)  = z"(xi)  + 0(hu(p+k)) 


Analogously  construct  operators  Rfa:  -*•  R such  that 

Rb(Ah  z)  = z 1 (b ) + 0(hu(p+k)) 


Def i ne 


28 


Vr(t) 


r = 1,  2,  ....  u - 
then 


/5  0 ’ “ 

- f(xi,  £.) 
\RbU)  * 6 


i=l,2 


^h,r(Ah  z)  = Ah{p(z)}  + 0(hM(p+k)) 

Hence  all  <j>h  r,  r = 1,  2,  v - 1 are  identical  and  they  are  consistent 
of  order  u(p  + k)  with  F.  From  the  formulas  it  is  obvious  that  we  do  not 
need  to  use  this  trick  for  the  approximations  to  z 1 (b ) but  can  use  operators 

Rb  r:  Eh  R such  that 


b ,r  h 


z)  = z'(b)  + 


M 

l 


v=( r+1 ) *p 


rv 


0(hM+1) 


h ,r‘ 


in  the  definition  of  <j> 
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3.  Approximation  of  linear  functionals 

To  construct  the  perturbation  operators  ^ and  ^ . , i = 1,  2,  ... 
in  the  applications  of  section  4 we  need  in  many  cases  formulas  that 
approximate  the  value  and  the  value  of  the  derivatives  of  a function  z for 
a given  argument  by  linear  combinations  of  the  components  of  Ah  z. 

Such  formulas  have  been  examined  in  Ballester,  Pereyra  (1967), 
Bjorck,  Pereyra  (1970),  Galimberti,  Pereyra  (1970),  (1971),  and  Pereyra 
(1973),  so  we  simply  refer  to  those  papers  and  introduce  some  notation 
that  will  be  useful  in  section  4. 

Let 

E = CS(a,  b)  , Eh  = RP+1 


Define  a^:  E -*■  Eh 

Ah  Z = [ Z ( Xg ) » z(  X-j  ) , ...,  Z ( Xp ) ] 

Xj  € [a,  b]  j = 0,  1 , . . . , P 


for  z € E. 

Define  the  linear  operators  (depending  on  h) 

x £ [a,  b] 

D^m:  Eh  - K v = 0,  1,  ....  P 

m=l,2,  ...,P+l-v 


according  to 


P 

l 

r=0 


a 


V 


C 


r 


(most  of  the  a (x)  may  be  zero)  such  that 
v,m,r 

5 

D^m(Ah  z)  = z(v)(x)  + Z g.(z)(x)  hj  + 0(hS+1) 

j=m  J 


for  z ^ E. 
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v m 

Note:  As  P is  finite  the  order  of  consistency  m of  the  operator  D ' with 

z^(x)  cannot  exceed  P + 1 - v. 

When  v,  m and  x are  fixed  the  constants  m r(x);  r = 0,  1,  ....  P can  be 
obtained  as  the  solution  of  a Vandermonde  system  of  linear  equations.  In 
Bjorck,  Pereyra  (1970)  an  efficient  algorithm  for  the  numerical  solution  of 
such  linear  systems  is  given. 

For  a Vandermonde  system  with  n unknowns  this  algorithm  requires 
approximately  3n^  operations.  To  get  D^,m  consistent  of  order  m with 
z^(x)  we  only  need  to  have  m + v of  the  coefficients  o / 0,  i.e. 

V jin , r 

the  Vandermonde  system  we  have  to  solve  has  only  m + v unknowns.  Further 
most  of  the  operators  D^’m  needed  have  x = x^  where  x..  is  a gridpoint. 

Once  we  have  found  the  operator  D^’m  for  one  gridpoint  we  can  easily  get 
the  operators  D^,m  for  most  of  the  other  gridpoints  (if  the  gridpoints  are 
equidistant)  without  solving  a system  of  linear  equations.  In  conclusion, 
the  amount  of  work  needed  to  find  these  operators  is  in  general  negligible 
compared  to  the  amount  of  work  involved  in  solving  the  operator  equations 

4>h ( n ) = 0 and  <}>h(nE)  + <f>E(n)  = 0 

or 

4>h(ni)  + I 4>h  i(ni  ^)  = 0 , i = 1,  2,  ... 

n v=l  n>1 

For  some  examples  of  operators  of  this  type  see  note  6 after 
theorems  1 and  2. 

For  functions  of  several  variables  the  formulas  above  can  be  used 
for  each  of  the  variables  or  one  can  use  more  compact  approximation  .formulas 
taking  linear  combinations  of  the  function  values  at  all  the  gridpoints. 


✓ > 
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To  distinguish  between  the  different  partial  derivatives  we  adjoin 
to  D the  variable  that  we  differentiate  with  respect  to,  e.g. 

v m ^ 

Dx’  is  consistent  of  order  m with  (— )(x,  y) 

,y  3xv 

v m „ v 

DTx]y,z,t  iS  consistent  of  order  m with  (-^-)(x,  y,  z,  t) 
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4.  Applications 

The  main  emphasis  of  these  applications  is  the  construction  of  the 
perturbation  operators  <t>^  and  4>h  ^ , i = 1 , 2,  . . . . The  operators  all 
correspond  to  well  known  discretizations  from  the  literature.  We  do  not 
claim  that  the  methods  chosen  are  the  best  possible  for  the  actual  problems; 
rather  we  have  tried  to  use  methods  that  are  fairly  well  known  and  in  some 
cases  widely  used. 

Although  the  examples  of  this  section  are  mainly  given  as 
illustrations  of  how  to  apply  the  gene  }1  ideas  of  section  2 to  different 
classes  of  problems,  we  believe  that  the  algorithms  for  iterative  improve- 
ment of  elliptic  partial  differential  equations  and  for  iterative  improve- 
ment of  integral  equations  compare  very  favorably  with  existing  algorithms 
for  these  classes  of  problems. 

The  questions  of  the  smoothness  of  the  expansions  of  the  global 
discretization  errors  for  the  basic  discretizations  have  not  been  examined 
in  detail.  For  the  different  applications,  wherever  possible,  reference 
to  known  results  on  such  expansions  are  given,  while  in  other  cases  we 
discuss  very  briefly  the  possible  existence  of  such  expansions.  These 
questions  need  further  study. 

In  the  theorems  of  section  2 we  make  some  assumptions  on  the 
perturbation  operators  and  ^ .. , i = 1 , 2,  . . . . The  validity  of  these 
assumptions  can  easily  be  checked  in  all  the  applications. 

In  future  studies  on  each  of  the  applications  of  interest  to  us 
a more  detailed  analysis  will  be  given. 

Some  numerical  results  are  given  in  most  of  the  subsections  below. 
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4.1  Initial  value  problems  for  ordinary  differential  equations 
Consider  the  scalar  differential  equation 
r = f(y)  ; y(0)  = a x 6 [o,  T]  f € CS(R) 

The  use  of  an  autonomous  problem  is  only  for  convenience  in  notation, 
as  is  the  restriction  to  a scalar  differential  equation.  All  the  results 
can  easily  be  generalized  to  non-autonomous  systems  of  ordinary  differential 
equations.  The  existence  of  smooth  asymptotic  error  expansions  for 
discretization  methods  for  initial  value  problems  for  ordinary  differential 
equations  is  discussed  in  Stetter  (1965),  (1973). 

In  our  operator  formalism  the  problem  is 


F: 


E 


E 


0 


E = r(0,  T) 


E°  = RxC(0,  T) 


max  h^(x) 


L|iE  OsxsT  v^Q  v! 

l9llEo  - 1^1  * o“*T  law 


for  g = 


1 


F(z)  = 


g(x)  0 s x * T 
z(0)  - a 

z'-f(z)  0 { x { T I 


Consider  Euler's  method 


y.+1  - y.  = h-f(y.)  1-0,1 N - 1 h = T/N 

This  simple  method  is  chosen  for  illustrative  purposes.  Although  the 
analysis  is  simple  for  this  method  the  techniques  used  to  construct  the 
perturbation  operators  are  applicable  to  more  complicated  methods  and  more 
complicated  problems. 
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Euler's  method  in  our  operator  formalism  reads: 

nN  + 1 


A.:  E - E. 

h h 


E = R 
h " 


llcll 

Ah  z (z(xq)*  z(x]),  ....  z(xN)) 


= maX  Ir 

Eh  OsuN  S' 


Let 


'2 


X 

II 

zr 

11 

O 

...  N ; 

A0  • F°  r° 

V E - Eh  ; 

E°  = Rx  RN 

1 91 

g = 

for  g = 

i‘0,1 N-l| 

9' 

g(x) 


.ill,  max  1 
r0  |Y  ' OsisN-1  *Yi 


and 


M) 


V Eh^Eh 


f«o  - a 


S+i  ' «i 


- f(5l) 


i =0,  1,  ....  N - 1 


Note  that 


z(xi  + l ) ‘ z(xi}  „ S-l  z(j+1)(x.) 


— - f(z(x.))  = z'(x.)  - f(z(x.))  + z -rrv  if-  hJ 

j=l  u l;- 


for  any  z £ E 
so 


S-l 


VAh  = Ah{F(z)  + E Mz)  hJ>  + 0(hb) 

j=l  J 


s x s T I 


+ 0(hS) 


I 


,(j+l) 


(I 


a 


LL 


0 < x s T 


Now  construct  4^  and  $h  v>  v = 1,  2,  ...  . We  will  give  some  different 
E 


operators  <t>h  v in  order  to  illustrate  some  useful  techniques  for 
construction  of  these  operators. 

1 . 


;U) 


vP't 

o 

p 

-e2  * «{,  - 

3^0 

2h 

Ci+1  " Ci-1 

\ « i 

- f 

- 


2h 


f(fM) 


i = 1,2 N - 1 


Note  that 

-z(x2)  + 4z(x1)  - 3z(xg) 


2h 


- f(z(xQ) ) = z* (xQ)  - f (z(xQ) ) 


S-l 


4 - 2j+1  _(j+l), 


* jS2  irrnTr2W  'J(xol'hJ  + 0(h  1 


and 


z(xi+])  - z(x._1) 


2h 


- f(z(xi))  = z'(x.)  - f(z(x.)) 


for  i =1,2,  ...,N-1 

Both  the  expansions  are  valid  for  any  z £ E.  Hence,  for  any  z £ E 

*h(Ah  z)  = Ah{F<z)}  + 0(h2) 

but  we  do  not  have 
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S-l 


♦hUh  z)  = a°{F(z)  + Y f.(z)  hj}  + 0(hS) 

j=2  J 


V h ” “h 1 

because  of  the  different  expansions  for  Xq  and  the  rest  of  the  gridpoints! 


2.  By  extending  the  operators  so  an  approximation  to  the  solution  at 
x_.|  = -h  is  computed  by  the  formula 


y0  - *-} 


,T^-  f(y.l)  -0 


we  can  define 


An  ' a 


*hu)  - 


Yl  ' ei-l 
2h 


f<«i> 


i = 0,  1,  ....  N - 1 


Alternatively  we  can  expand  the  operators  so  that  an  approximation  at 
XN+1  = xn  + h is  comPuted  by  the  formula 

yN+l  " yN 


f(yN)  = 0 


and  define 


;(c)  = 


40  ' “ 


ii±2  - 4?i+l  ' 3?i 


2h 


In  both  cases  we  have,  for  any  z £ E 


f<Y 


i = 0 , 1 , 


N - 1 


.E 


.0, 


S-l 


Vh(Ah  z)  = a”{F(z)  + Z f.(z)  hJ)  + 0(hb) 

The  operators  f.  will  of  course  be  different  for  the  two  different  operators 
j 

By  extending  the  numerical  solution  even  further  outside  the  interval 

of  interest  in  this  way  one  can  easily  construct  operators  , 

n » v 
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v = 1 , 2,  ...  such  that 

n S-l  , s 

4>h  v(Ah  z>  = Ah{F(z>  + 1 fv  i(z>  h } + 0(h  ^ 

h,v  h h j.(v+l)  V’J 

For  these  extended  operators  the  spaces  E,  E°,  Eh  and  E^  must  be  modified 
in  order  that  the  theorems  of  section  2 be  applicable.  Slight  modifications 
of  the  theorems  may  also  be  necessary. 


M)  = 


- a 


- fUi) 


i = 0,  1,  ...,  N - 1 


v = 1,  2,  ... 


Here  D^,S  are  operators  of  the  type  discussed  in  section  3.  For  any  z € E 


t»hjV(Ah  z)  = A°{ F( z)  1 + 0(hS) 


v = 1,  2,  ... 


C0  - ^ 


;(e)  = 


*i+l  " Ci  1 


\ — • i [f^i+i>  + 

Note  that  for  y £ E such  that  y'  = f(y) 


i =0,  1 N - 1 


y(xi+1)  - y(x.) 


- \ Cf(y(xi))  + f(y(xi+1 ))] 


S y^v^(xi)  t r 


i z f'r;(y(xi))(  z — rr^-b  > 
1 r=l  1 t=l  z' 


S y(v)(x.)  u , , S ,, 

= y'(x.)  - f(y(x  ))  + z — — j— h 1 - - r ‘ 

v=2 

+ 0(hS) 

S-l  i s 

= y'fxJ  - f(y(x.))  + z g,(y(x.))*hJ  + 0(hb) 

1 1 j=2  3 1 

The  absence  of  a term  g.|(y(x..))h  is  due  to  the  fact  that  for  the  elements 
y that  we  consider  we  have  y"  = f'(y)  y1! 

Note  that  the  expansion  above  is  not  valid  for  arbitrary  z £ E, 
but  only  for  z € E such  that  z1  = f(z).  In  this  case  we  have  to  rely  on 
theorem  2 while  for  the  previous  perturbation  operators  we  can  rely  on 
both  theorem  1 and  theorem  2. 

Note  that  for  y £ E 3 y'  = f(y), 

S-l 

*jj(Ah  y)  = A°{F(y)  + z f<(y)  hj>  + o(hs) 


j=2 


Further  note  that  any  method  of  second  order  or  more  can  be  used  to  construct 
the  perturbation  operator  ^ . The  main  advantages  with  our  perturbation 
operator  (which  is  based  on  the  trapezoidal  method)  is  that 

1)  no  new  values  of  f(y)  need  to  be  computed  to  get  the  perturbation 

2)  the  perturbed  solution  at  x..  can  be  computed  directly  after 
the  calculation  of  the  unperturbed  solution  at  xi 

3)  the  basic  discretization  formulas  for  and  are  both 
one-step  formulas. 

For  general  linear  multistep  methods 
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k 


I 

v=0 


av  Vv 


- h 


k 

z 

v=0 


*v  Vv  = 0 


one  can  e.g.  use  the  following  basic  discretization  formula  for  the 
perturbation  operator 


k 


z 

v=0 


rn+v 


- h 


k 

Z 

v=0 


n+v 


(We  assume  that  the  global  discretization  error  for  the  solution  obtained 
by  the  linear  multistep  method  has  a smooth  error  expansion;  this  may  not 
be  true  in  many  cases!)  Remember  that  the  maximal  order  of  a stable  linear 
multistep  method  of  stepnumber  k is  k + 1 when  k is  odd,  and  k + 2 when  k 
is  even  (see  e.g.  Lambert  (1973)).  By  choosing  a*,  B*,  v = 0,  1,  ....  k 
judiciously  one  can  get 


k 

z 

v=0 


'(Vv>  - 


k 

z 

v=0 


6*  f(y(xn+v))  = y'(xn)  - f(y(xn))  + 0(h2k) 


for  y£  E such  that  y'  = f(y). 

One  can  of  course  also  use  perturbation  operators  of  the  type  discussed  in 
point  3 above.  The  basic  discretization  operator  for  here,  however,  has. 
the  same  stepnumber  k as  the  basic  discretization  operator  for  <j>h  while  the 
stepnumber  for  the  basic  discretization  operators  in  point  3 may  be 
considerably  larger. 

No  study  of  the  practical  implementation  of  these  ideas  for  initial 
value  problems  for  ordinary  differential  equations  has  been  undertaken  yet. 


5.  As  a further  illustration  we  choose  a different  operator 
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& Eh’E°  ; 

eS- 

R x 

rn 

1 g1 

/g1 

\ 

"h  3 ’ 

for  g 

= 

y9(xi> 

i ■ 1,  2,  ....  N I 

^g(x) 

0 s x $ T 

Mill  0 ■ It’I 

max 

nun 

hi 

1 

The  operator  4>h  is  the  same  as  above,  i.e. 


*„U)  = 


'■o  - “ 


Si*l  ' Ei  */.  i 1=0,1 N - 1 


• f(5i> 

Note  that  for  i = 0 , 1 , . . . , N - 1 and  any  z £ E 

Z(Xi+1}h  - f(z(x.))  - z'(x.+1)  - f(z(x.+1)) 


s-1 


+ i 9i(z(x.+1 ) ) hJ  + 0(hs ) 
j=l  J 11 


Choose 


+hu)  ■ 


Eo  ■ “ 


- s 


i+1  ‘•i-1 

2h 


- f(^)  i = 1,  2,  ....  N 


then 


s-1 


*h(&h  z)  = VF(z)  + 1 fj(z)  h>3}  + °(hS) 

J 2 


J 
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Note  1 Some  of  the  results  of  this  section  carries  ov*r  to  the  case  when 
x.,  i = 0,  1,  N - 1 are  not  equidistant.  E.g.,  the  perturbations  of 
point  3 and  point  4 are  useful  in  such  cases. 

Note  2 Consider  again  the  scalar  differential  equation 

r - f(y)  = o y(0)  -a  ; x ( [0,  T]  ; f € CS(R) 

Note  that  this  problem  is  equivalent  to  the  Volterra  integral  equation 

x 

y(x)  - y (0)  - f f(y(s))  ds  = 0 
0 

Euler’s  method  for  the  original  differential  equation  can  be  viewed  as  an 
approximation  to  the  integral  equation,  namely 

yi+1  = + h f(y.)  = yi_1  + h f(yi_1)  + h f(y.)  = ... 


= a + h l f (y  . ) i =0,1,2,  ... 
j=0  J 

With  appropriate  definitions  of  the  spaces  E,  E®,  Eh>  E^,  the  operators 
A^ , A®  and  the  corresponding  norms  we  have 

x 

F(z)  = (z  - a - / f(z(s))  ds  , 0 5 X $ T) 

0 


*htt)  ■ 


a 


a 


1-1 

2 f(cJ 
j=0  J 


For  the  actual  computations  we  would  of  course  use  the  recursion  formula 

= « 

€i  = + h f(c1_1)  i = 1,2 N 


Note  that  is  consistent  of  order  1 with  F,  i.e.  with  the  integral 
equation. 
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Define 


♦h<«>  ■ 


^0  ' a 


- h T.  a..  f(c  ) 

j=0  1J  J 


1 = 1,  2,  ....  N 


where 


ai0  = aii  = 1/2 


a>.jj  = 1 j = 1 , 2,  . . . , i - 1 


For  the  actual  computations  we  would  use  the  formula 

[^(£)]0  -£0-<. 

= C*h(i:):1i-1  * '•i  - £j-l  - jr  * f(£,.,>]  1 ■ 1,  2,  3 

Note  that  is  consistent  of  order  2 with  F.  Other  perturbation  operators 
could  be  used! 

The  advantage  of  this  point  of  view  is  that  for  the  maximum  norms 


ill  . WdX  | 

r0  ^ 1 OsisN  *Yi 


where 


and 


where 


Y = 


Yi 


i = 1 (1 )N 


€ E 


IUIlEh  = 0^N  Ui 


C = (^  , i = 0(1)N)  € 


we  have 
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I lye1)  - ys2)ii  0 * c*  i u1 . s2n, 

Eh  h 

II^U1)  - *[;U2)1I  0 * ce‘i|c1  - ^llE 

Eh  h 

i.e.  k = 0 in  condition  e)  of  theorem  1,  while  for  Euler's  method  directly 
we  have  k = 1 for  the  maximum  norm.  The  same  result,  i.e.  k = 0,  may  be 
obtained  for  Euler's  method  directly  if  we  use  Spijker's  norm  for  E^  (see 
Stetter  (1973),  section  2.2.4,  p.  81-84)  which  in  our  case  reads 


for 

i » 1(1  )N  - 1 

Note  the  relation  between  Spijker's  norm  for  E^  corresponding  to  Euler's 
method  and  the  maximum  norm  for  E^  corresponding  to  the  integral  equation 
formulation. 

The  reformulation  of  the  original  problem  as  an  integro-di fferential 
equation  may  be  a useful  trick  for  other  methods  for  initial  value  problems 
for  ODEs  and  initial  boundary  value  problems  for  PDEs. 

The  ideas  presented  in  this  note  need  further  investigations  and 


- SS* 


are  included  mainly  as  an  illustration  of  a way  to  increase  the  range  of 
problems  and  methods  for  which  the  theorems  may  be  useful. 
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^ ^ Two-point  boundary  value  problems  for  ordinary  differential  equations 
Consider 

y"  = f(x,  y,  y 1 ) a-e*x*b  + e 

r^yfa),  y'(a))  = 0 r2(y(b),  y * (b) ) = 0 

Assume  that  the  functions  f,  r^  and  r^  are  such  that  the  boundary  value 
problem  has  a unique  solution  which  is  2s  times  differentiable. 

Define 

F:  E - E° 


E = C2s(a 


» b + e) 


II  . max  Is  |z(,)(x)| 
"E  * v=0  v! 


E =RxC(a,b)xR 
_o  = It1!  + |y2|  + mxX  |g(x)| 


for 


Introduce 


1 

Y 

g(x) 

2 


a-esx$b  + 


F(z)  = 


\ 

I r^zU),  z'  (b) ) 

z"  - f(x,  z,  z’) 
^ r2(z(b),  z ■ (b ) ) 


€ El 


a-esx$b  + 


V E^Eh 


E = RN+2 

n " 


Ah  2 ” £z(xq)»  •••»  z(x^+-|)] 


where  x^ 


a - j + 1-h,  i = 0.  1,  ....  N + 1;  h * 


. 0 r0  r0 
V E -*Eh 


E°  = R x RN  x R 


Defi ne 


V Eh-Eh 


. ,S  + *0  S - «o. 

ri(— T~  * — h — } 


S'+l  ' 2 S'  + S-l 

u2 


f(x  r " S-K 

1 i ’ S'  • 2h  J 


i = 1,  2,  ....  N 


\ ..  Sn+i  + S S+i  ‘ Sk 

\ r2  2 * h > 

Assume  that  f,  and  are  such  that  t)>h  is  stable.  is  consistent  of 
order  2 wi th  F. 

The  existence  of  smooth  error  expansion  is  discussed  in  Stetter 
(1965),  Pereyra  (1968). 


Construction  of  perturbation  operators 


1 . Direct  approach . 


.E.  F F0 

V Eh"Eh 


rl(Da’4(c)  * °a’4^)) 


^(C)  = o;*  (0  - f(xif  cr  d;;Sc)) 

\r2(Dg’4(c),  dJ‘4(C)) 


i = 1,  2,  ...,  N 


where  D^’  : E^  -*•  R are  operators  of  the  kind  discussed  in  section  3. 


Define 
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3.  Extending  the  solution  outside  the  interval  [a  - e,  b + e]  (from  an 
unpublished  paper  by  H.  Keller,  V.  Pereyra,  communicated  by  V.  Pereyra). 

To  be  able  to  use  symmetric  and  identical  formulas  to  approximate  z^(x^), 
v = 1,  2 at  all  gridpoints  xi  (including  the  gridpoints  close  to  the  boundary) 
one  can  extend  the  numerical  solution  to  "artificial"  gridpoints  outside 
the  interval  (a,  b)  by  the  basic  formula 
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Vl  ' 2*n  + Vl 


- f(x  , y ) = 0 
' n’  •7rr 


Sufficiently  many  exterior  points  are  introduced,  so  symmetrical  and  identical 

approximation  formulas  for  the  first  and  second  derivative  of  any  z € E can 

be  used  for  all  gridpoints  needed.  Now  introduce  <p.  . , i =1,2, 

n , i 

Vl(e)  ■ | - f(Xj.  Cj.  j ■ 1.  2 N 

\ '-2(''2'2<’*,)(0,  o!'2(,t|)(5))  / 


j(Ah  z)  = A°{F(z)  + r f..(z)  hJ}  + 0(h2s_1) 
j=(i+l).p  ,J 


with  a proper  definition  of  A,  and  A?. 

h h 

Numerical  results 

The  following  special  cases  of  (1) 

W y"  - f(x,  y,  y')  = 0 

y(a)  = a y(b)  = e 
(B)  y"  - f(x,  y,  y')  = 0 

y(a)  = a y 1 (b)  = 6 

were  discretized  by  formulas  analogous  to  those  above.  As  the  boundary 
conditions  for  these  cases  are  simpler  than  in  the  general  case,  somewhat 
different  discretizations  are  used  here: 


3 a + i-h  i = 0,  1,  ....  N 


b - a 


- « = ° 
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’■i  + 1 " * Ci~1  . f(x  , Ci+1  ~ gi-K  n 

u 2 (xi  ’ ci  ’ 2h  ^ ■ 0 


i = 1,  2,  ....  N - 1 


«N  ' 6 = 0 


xi  = a + i-h  1 * 0,  1,  ....  N h = 


2N  - 1 


CQ  - a = Q 


^ * ;j-'  - fix  5 hn  : 

h2  ' Xi  ’ 5i’  2h 


) = 0 


CN  " CN-1 


6 = 0 


The  perturba on  operators  <f>^  we  use  are 


co  • a 


*h(c)  = Dx!4(c)  ' f(xi’  V Dx!4^))  1 = 1,  2,  N - 1 


*N  - 6 


*0  - a 


*h(c)  = V^)  ‘ f(V  V Dx!4(c))  1 = ]>  2«  •••»  N - 1 


\dJ’4(0-6  J 

The  systems  of  non-linear  equations  obtained  by  the  discretization 
were  solved  by  Newton  iteration  and  the  initial  approximation  0 was  used. 
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The  iterations  were  terminated  when  the  last  correction  in  the  iterations 
was  less  than  e. 

The  perturbed  problem  was  solved  by  a modified  Newton  iteration, 
where  the  LR-factorization  of  the  final  iteration  matrix  obtained  in  the 
solution  of  the  unperturbed  problem  was  used.  As  initial  approximation 
the  solution  of  the  unperturbed  problem  was  used  and  the  iterations  were 
carried  on  until  an  estimate  of  the  iteration  error  was  either  less  than 
one  tenth  of  the  estimated  discretization  error  or  less  than  e.  All 
quantities  were  measured  in  the  maximum  norm. 

In  the  table  below  the  entries  of  the  column  "Number  of  iterations" 
stand  for  the  number  of  iterations  for  the  unperturbed  problem/number  of 
iterations  for  the  perturbed  problem. 

The  following  differential  equations  were  solved: 

(Al)  y"  = \ (y  + x + l)3  y(0)  = y(l)  = 0 

exact  solution:  y(x)  = 2/(2  - x)  - x - 1 

(Ciarlet,  Schultz  and  Varga  (1967)) 

(A2 ) y"  = y3  - sin(x)  * (1  + [sin(x)]2)  y(0)  = y(ir)  = 0 

exact  solution:  y(x)  = sin(x) 

(Pereyra  (1968)) 

(A3)  y"  = - 1 - 0.49(y’)2  y(0)  » y(l)  = 0 

exact  solution:  y(x)  = ^ ln( Vcos\^r2Sy  ") 

(Bellman,  Kalaba  (1965)) 

(Bl)  y"  = \ (y  + X + l)3  y ( 0 ) = 0 y'(l)  = 1 

exact  solution:  y(x)  = 2/(2  - x)  - x - 1 
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The  problems  were  solved  in  single  precision  on  an  IBM  360/75  with 
N = 10  and  e = 10"6. 


1 

i Problem 

| |Actual  error| | 

| |Estimated  error| | 

Number  of 
iterations 

mo 

A1 

7.3  • 10"4 

1 

O 

4/2 

1.7  • 10'1 

A2 

2.38  • 10'3 

2.36  • 10"3 

6/2 

i 

A3 

1.059  • 10'4 

1.060  • 10'4 

3/2 

1.3  • 10'1 

B1 

1.12  • 10"3 

1.04  • 10"  3 

4/2 

o 

1 

Obviously  one  can  get  very  good  error  estimates  at  a low  cost  for  these 
problems.  All  the  problems  are  very  simple  with  fairly  rapid  convergence 
of  the  Newton  iterations  for  the  original  discretization.  In  more 
difficult  cases  where  it  may  be  essential  to  have  a good  initial  guess  for 
the  solution  the  quotient  of  the  amount  of  work  for  the  unperturbed 
problem  to  the  amount  of  work  for  the  perturbed  problem  may  be  much 
bigger. 

No  experiments  have  been  done  with  iterative  improvement  for  this 
class  of  problems. 

4.3  Two-dimensional  elliptic  boundary  value  problems 

We  will  discuss  two  classes  of  problems  on  rectangular  regions. 

All  problems  discussed  are  assumed  to  have  unique  smooth  solutions, 
however  some  numerical  results  for  a problem  with  a non-smooth  solution 
will  be  presented. 

The  existence  of  smooth  error  expansions  for  discretizations  of 

elliptic  problems  are  discussed  in  Volkov  (1957),  Stetter  (1965),  Hofman 

\ 

(1967)  and  Pereyra  (1970). 
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4.3.1  Problems  non-linear  in  u only 
Consider 

p 

V u = f ( x , y,  u) 


u = h(x , y) 


(x,  y ) € n 
(x,  y)  £ an 


where 


v2  = iL  + si- 

ax2  ay2 

a 10  $ x 5 li  0 s y j 1} 

an  the  boundary  of  n 

Or 

Assume  that  f and  h are  such  that  the  solution  u 6 C (n). 
In  our  operator  formalism  we  have 

.0 


E = C2s(n) 


F:  E -+  E 


-0 


E = C(n)  x C(an) 


with  the  norms 


max 


Is  J.  l iiLzhj-jli 

1 v=0  v!  u=0  axu  a yv'w 


where 


r0  = (x^y)  € n |gl(x’  y)|  + (x,"y)  £ |g°(x’  y)| 


g1 (x,  y)  (x,  y)  £ a 

g°(x,  y)  (x,  y)  £ an 


F(z)  = 


v*  z - f(x,  y,  z) 
z - h(x,  y) 


(x,  y)  £ n 
(x,  y)  € an 
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Introduce 


and 


V 


E.  = R<"+')2 


Ah  2 = (z<xi’  yj) 

x.  = i*h 


i = 0(1 )N > j = O(l)N) 
i = 0( 1 )N 


Yj  = j-h  j = 0(1 )N 
h = 1/N 

I U I lE  = max  UiJ 

h OsisN  1J 
Osj^N 


A 


0 

h 


g 


for  g € E° 


as  above. 


F°  - r(n-1  )' 
Lh  K 


x RN+1  x RN+1  x RN  x RN 


9°(x0,  y.) 
g°(xN,  yj) 
g°(xi,  yQ) 

g0(xi’  yN> 


ri  = 1(1 )N  - 1 
{j  = 1(1)N  - 1 

j = 0(1 )N 
j = 0(1  )N 
i = 1 (1  )N  - 1 
i = 1 (1  )N  - 1 


0 = tnax[  | y!  . | ;i  = 1(1)N  - 1,  j = 1 (1  )N  - 1, 
h 


= 0(1  )N,  |Yf|;j  = 0(1  )N, 

|y?3|;i  = 1(1)N  - 1,  |y°4|;1  = 1 (1  )N  - 1] 


where 
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Defi ne 


*h U)  = 


. 1 

i = 1 (1  )N  - 

Yij 

j = Kl  )N  - 

01 

Yj 

j = 0(1  )N 

Y = 

02 

YJ 

j = 1(1  )N 

03 

Yi 

1 = 1(1  )N  - 

04 

_Yi 

i = 1(1  )N  - 

V Eh-Eh 

?i+lj  " ci 

j+1  ' 

4ti.i + H-1.1 + E 

€ E 


1JZ1  ' f(x1’  yj’  Cij} 


i = 1 ( 1 ) N 

C0j  ' h^V  yj^ 

j = 0(1)N 

CNj  “ h^V  V 

j = 0(1 )N 

eiO  ” h^xi * y0^ 

i = 1 (1  )N 

CiN  " h^xi * yN^ 

i = 1 ( 1 ) N 

One  can  easily  verify  that 


2s-2 


♦hUh  z)  = a[J{F(z)  + I f.(z)  hj}  + 0(h2s_1) 


j=2 


Let 


♦h(«)  = 


,(jr) 

Hj' 

c0j  - h(xQ,  y j ) 
SNj  - h(xN,  y.) 

eiO  " h^xi  ’ yO^ 
CiN  - h(x.,  yN) 


i = 1(1 )N  - 1 , j = 1(1 )N  - 1 
j = 0(1  )N 
j = 0(1 )N 
i = 1 ( 1 )N  - 1 
i = 1(1 )N  - 1 
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where  <|».  . : E.  ->•  R. 

l j h 

We  will  consider  three  different  alternatives  for 

' J 


1 . 


where  DX  and  DY  are  operators  of  the  kind  discussed  in  section  3.  It  is  easy 
to  verify  that 


^(Ah  z)  = A°{F(z)}  + 0(h4) 


for  all  z £ E . 


2-  ♦ijtO  • (ei.i  j-i  + ^i-ij  + «i.i  j+i  + *«iM  - 20^ij  + 4^ij+i 
+ Ci+1  j-1  + 4ci+l j + C1+l  j+l)/(6h2) 

- tJ  [f(V!  * yj-i  * ?i-l  .j-1  ^ + 4f(xi-l  ’ yj’  ci-l  ,j^ 

+ f(xi-r  yj+r  €i-i.j+i)  + 4f(V  yj-T  «u-i) 

- 20f(xr  y j » 5lj)  + 4f(x.,  y.+1,  CiJ+1) 

+ fUi+1.  yj_i » ci+i ,j-i ^ + 4f(xi+r  yj#  ci+i,j) 

+ f(xi+r  yj+r  ci+i,j+i)]  ‘ f(V  yr  cij) 

One  can  show  that  for  the  solution  u of  F(u)  = 0 

- n2 


0ij(Ah  u)  = u(xi,  yj)  - f(xi,  y. , u(x.,  y j ) ) + 0(h  ) 


so 


4>|jUh  = + 

for  the  solution  u of  F(u)  = 0.  In  this  case  we  have  to  rely  on  theorem  2, 
while  for  perturbations  1 and  3 we  can  rely  on  both  theorem  1 and  theorem  2. 
Further  one  can  show  that 

2s-2 

■*(a.  u)  = a°{F(u)  + z f.  ,(u).hj}  + 0(h2s_1) 


V“h 


j=4 


E ,j 
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2 

We  have  used  the  ninepoint  operator  vg  to  approximate  the  Laplace 
2 

operator  v . For  notation  and  the  results  below  see  e.g.  p.  321  in  Bjorck, 
Dahlquist  (1973). 

We  have 

vg  u = v2  u + - — + °(h4) 

But  v u = f(x,  y,  u)  so 

v2  u = f(x,  y,  u)  + h2  v2  f(x,  y,  u)/12  + 0(h4) 

p 

Thus  if  we  use  the  ninepoint  operator  to  approximate  v f(x,  y,  u)  we  get 

2 

v2  u - f(x,  y,  u)  - v9  f(x»  y»  u)  = °(h4) 


which  gives  the  result  above, 
then 


rj 


rj 


^(Ah  z)  = A°{F(z)}  + 0(hq) 


so  if  q * 2*(v  + 1 ) we  can  take 

m3  x 


9h,v  = % 

and  make  v iterative  improvements, 
max 


v = 1 , 2 , ....  v, 


max 


Numerical  results 

The  problem  above  with  f = 0 and  h(x,  y)  = sin(x)  sin  h(y) 

2 2 

+ cos  h(x)  cos  y - (x  - y )/2  was  discretized  as  above  (Kronsjo,  Dahlquist 
(1972)).  The  system  of  linear  equations  obtained  was  sclved  iteratively 
with  SOR  with  optimal  relaxation  parameter,  and  to  get  an  initial  approxi- 
mation we  interpolated  the  boundary  values  linearly.  The  iterations  were 
terminated  when  the  residual  was  less  than  t. 
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The  perturbation  operators  of  both  1 and  2 above  were  used.  The 
perturbed  problem  was  solved  with  the  same  method  but  now  we  used  the 
solution  from  the  unperturbed  problem  as  initial  approximation  and  the 
iterations  were  terminated  if  an  estimate  of  the  iteration  error  was 
either  less  than  y times  the  estimated  discretization  error  or  less  than  e. 
The  errors  were  measured  in  the  maximum  norm. 

The  calculations  were  performed  in  double  precision  on  an  IBM  360/75 
with  N = 10  and  some  different  values  of  e and  y.  The  entries  in  the  column 
"Number  of  iterations"  stand  for  the  number  of  iterations  for  the  unperturbed 
problem/the  number  of  iterations  for  the  perturbed  problem. 


Perturbation 

c 

Y 

Actual 

Error 

Estimated 

Error 

Number  of 
Iterations 

1 

0., 

1.5M0*4 

1.40-10'4 

32/7 

10'10 

1.5M0'4 

1 .50*  TO"4 

44/21 

10'4 

10'4 

1 .27*1 0-4 

5 .3-1 0-5 

21/1 

o 

1 

cn 

10‘4 

1 .47* 10*4 

1 . 36- 1 0~4 

24/8 

IQ'6 

*3* 

1 

o 

1.5M0'4 

1 . 50 - 1 0~4 

29/14 

2 

10-7 

10‘4 

1 .51 • 10"4 

1 .50- 10-4 

32/18 

10"7 

0.1 

1.5MQ'4 

1 .40- 10-4 

32/7 

10-10 

*3- 

1 

o 

1 .5M0'4 

1 .50- 10"4 

44/21 

Note  that  there  is  no  difference  between  the  results  for  the  two  different 
perturbations.  Further  note  that  to  get  a reliable  error  estimate  we  must 
get  the  iteration  error  sufficiently  small.  The  estimation  algorithm  can 
only  estimate  the  discretization  error,  not  the  iteration  error.  In  fact, 
there  is  a danger  of  amplification  of  the  iteration  errors  when  we  apply 
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the  operator  <j>^  to  the  solution  of  the  unperturbed  problem.  Note  also  that 

we  can  save  some  work  by  terminating  the  iterations  for  the  perturbed  problem 

early,  i.e.  by  choosing  a larger  value  for  y. 

Analogous  results  were  obtained  for  Laplace  equation  with  derivative 

boundary  conditions  on  some  of  the  sides  of  the  unit  square. 

The  same  problem  was  solved  with  iterative  improvement,  using  the 

perturbation  operators  $,  , v = 1,  2,  ...  according  to  point  3 above.  We 

n , v 

-14  -9 

used  e = 10  , y = 10  and  some  different  values  of  N and  q.  The  tables 

below  give  the  maximal  errors  in  the  successive  iterates  for  the  iterative 
improvement. 


N = 10 
q = 8 

Iteration  number 

0 

1 

2 

3 

4 

5 

6 

Error 

1.5  10"4 

7.8  10'9 

1.4  10‘9 

2.7  10’10 

5.6  10"11 

Number  of 
Iterations 
in  SOR 



52 

40 

37 

31 

26 

22 

N = 10 
q = 4 

Iteration  number 

0 

i 

2 

3 

4 

5 

6 

Error 

1.5  10'4 

1.2  10"6 

6.6  10"8 

6.9  10"8 

6.9  10"8 

6.9  10'8 

6.9  10-8 

Number  of 
Iterations 
in  SOR 



52 

40 

37 

30 

24 

20 



19 

N = 20 
q = 12 

Iteration  number 

1 

2 

3 

4 

5 

6 

Error 

3.8  10-5 

7.6  10‘8 

4.2  10'9 

9.0  10‘10 

2.7  10"10 

9.9  10'11 

4.0  10*11 

Number  of 
Iterations 
in  SOR 

102 

80 

62 

51 

48 

45 

43 
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N = 20 
q = 8 

Iteration  number 

0 

1 

2 

3 

4 

5 

6 

Error 

3.8  10"5 

7.6  10'8 

3.4  10  ’ 

5.2  10'10 

9.8  10'11 

2.1  10'11 

-12 

4.4  10 

Number  of 
Iterations 
in  SOR 



102 

80 



62 

51 

47 

43 

41 

Note  1 In  all  cases  but  one  (q  = 12)  we  have  made  more  iterative  improvements 
than  suggested  by  theorems  3 and  4.  In  the  theorems  it  is  implied  that  q 
should  be  chosen  such  that  q > (maximal  number  of  iterations  + l)-2.  Further 
note  that  we  improve  the  accuracy  of  our  solutions  even  after  the  suggested 
maximal  number  of  iterations  has  been  exceeded.  For  q = 4 e.g.,  only  one 
improvement  is  reasonable  according  to  the  theorems,  but  the  second  iterate 
has  a considerably  much  smaller  error  than  the  first.  This  astonishing  result 
is  a lucky  coincidence  in  the  construction  of  the  perturbation  operators 
, v = 1,  2,  ...  . For  almost  all  the  gridpoints  (all  but  the  points 

n , V 

closest  to  the  boundary  points)  the  formulas  we  use  are  of  order  6 (rather 

2 

than  of  order  4)  for  all  functions  u such  that  v u = 0.  A similar  result 
holds  for  q = 8.  If  <j>^  v(a^  u)  = 0(hc')  we  cannot  increase  the  order  of 
accuracy  of  our  solutions  by  making  extra  iterations,  however  the  "error 
constant"  can  in  some  cases  be  decreased  in  this  way. 

Note  2 In  the  cases  q = 8 and  q = 12  the  maximal  errors  in  the  last  three 

iterations  occur  close  to  the  boundary.  At  the  interior  points  the  errors 

are  much  smaller.  The  explanation  for  this  is  that  for  the  points  close  to 

the  boundary  we  have  to  use  non-centered  approximations  to  the  derivatives, 

while  at  the  interior  points  we  can  use  symmetric  approximations. 

Note  3 There  is  a certain  amplification  of  the  iteration  errors  when  *. 

n , v 

is  applied  to  the  solution.  The  larger  q is  the  larger  is  this 
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amplification.  Further  the  larger  q is  the  larger  are  the  "error  constants" 
for  the  formulas  employed  close  to  the  boundaries.  These  facts  may  explain 
why  we  get  better  results  wi th  q = 8 than  with  q = 12  for  N = 20. 


4.3.2  The  minimal  surface  equation 

As  an  exercise  in  how  to  proceed  in  a more  difficult  case  we 

consider 


f-  r -1  — iM.]  + !_  [ 

3x  L o o 9xJ  3y  L 

''*4*  “y 


1 |£]  --  o (X,  y)  e n 


/■.  2 2 

/I  + u + u 
x y 


u = (cos  h2(y)  - x2)1/2  (x,  y)(  3D 


where 


8 ! (0  { n 1,  0 j y J 1) 

In  vector  operator  notation  the  equation  reads 

V- [y ( I vu I 2)vu]  = 0 

where 


y(|vu|2)  = (1  + |vu|2)'1/2 
Note  that  the  exact  solution  is 

u(x,  y)  = (cos  h2(y)  - x2)1^2 

and  that  has  a singularity  at  x = 1,  y = 0.  The  problem  is  discretized 
and  solved  according  to  Concus  (1967)  in  the  following  way: 

Introduce  the  gridpoints 

x^  = i - h i = 0(1 )N 

yj  = j-h  j = 0(1  )N 

where  h = 1/N. 
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Approximate  the  differential  equation  by 

fij  = ‘ * 'rT*i  J<2uij  - ui+l.j  • 

* T,j*l(2ulj  - - “i .j+i1  * yT*i,jti(2uij  - Vi.j  ' 

= 0 , i = 1 (1  )N  - 1,  j = 1 (1  )N  - 1 

2 

where  yjj  = Y(|vu|yy)  denotes  y for  the  mesh  cell  with  center  (i  - 1/2, 

1 J 1 J 

j - 1/2),  which  is  evaluated  by  use  of 

* " u1-l ,j-l )2  * ' u1-l,j-l)2]- 

This  approximation  is  consistent  of  order  2 with  the  differential  equation. 
The  boundary  conditions  are  represented  by 

uij  = [cos  h2(yj)  - xj)1/2  (x. , y^)  € 3 n 

The  system  of  nonlinear  equations  above  are  solved  iteratively  by  computing 
k+1 

u^  , the  (k  + 1 ) th  approximation  to  u^  from 

f r k+1  k+1  k k -i 

uk+i  « uk . . u ij}un  »-•••»  UJ.-J y.y  un,n-i-< 

,J  1J  ' f i j r k+1  k+1  k k 7 ’ 

3U~  tull  ui  -1 , j ’ ui  j ’ •••’  UN  ,N-1  ^ 

where  u>  is  the  relaxation  parameter.  The  initial  approximation  was  obtained 
by  linear  interpolation  of  the  boundary  values  and  the  iterations  were 
terminated  when  the  last  correction  was  less  than  e. 

The  following  basic  discretization  formula  was  used  for  the 


perturbation  operator  4 
Introduce  the  notation 


£ ■ <£,j 


1 ■ 0(1)11.  j = 0(1)N) 


and 
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px(0  = 


py(c)  = 


xi’yi 


/I  + (DX1,4  U))2  + (dyJ’4  U))2 
Xi ,yj  Xi  ,yj 


DY1,4  (?) 

x . ,y . 


: /I  + (Dx'*\  (s)T  + (DY'*\  U))‘ 
Xi,yj  Xi,yj 


i = 0(1 )N 
j = 0(1 )N 


i = 0(1  )N 
j = 0(1  )N 


Here  DX  and  DY  are  operators  of  the  type  described  in  section  3.  Further 

* (r)  = DxJ’4  (px(O)  + DyJ’4  (py(e))  i = KDN  - 1,  j = 1(1)N  - 1 

1 j * i *y  j *y  i 


xi-yj 


Note  that  with  proper  definition  of  , 

z)  ■ C^r  ( 4 S It)  * = ( ■=--■>  * 0(h4) 


i j h 


— r—  £“]  + n ( * v ) + 

Vi  * u2  * u2  sx  » /i  + u2  - u2  »rJl  1>yJ 

x y X y 


is  the  basic  discretization  operator  for  the  perturbations.  The 
perturbed  problem  was  solved  with  the  same  iterative  method  as  the  unperturbed 
problem,  but  now  we  used  the  solution  of  the  unperturbed  problem  as  initial 
approximation  and  terminated  the  iterations  when  an  estimate  of  the  iteration 
error  was  either  less  than  y times  the  estimated  discretization  error  or 
less  than  t.  All  quantities  were  measured  in  the  maximum  norm. 

The  problem  was  solved  in  double  precision  on  an  IBM  360/75  and  we 
-8  -4 

used  e = 10  , y = 10  and  some  different  values  of  N and  w.  The  entries 

in  the  column,  "Number  of  Iterations,"  refer  to  the  number  of  iterations 
for  unperturbed  problem/number  of  iterations  for  perturbed  problem. 


N 

Actual  Error 

Estimated  Error 

Number  of  Iterations 

10 

1.7 

4.3*1 O-3 

8.9-10'3 

55/28 

40 

1.87 

2.3-10'3 

5.0-10'3 

147/71 
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The  maximum  error  was  obtained  in  the  vicinity  of  x = 1,  y = 0.  Due  to  the 
singularity  of  ux  at  that  point  we  cannot  expect  to  get  very  good  results 
(or  good  error  estimates)  in  that  region.  However,  for  the  rest  of  the 
unit  square  much  better  results  and  error  estimates  were  obtained.  Below 
the  results  at  some  representive  points  are  tabulated  for  N = 10  and 
u = 1.7. 


X 

* 

Solution 

Actual 

Error 

Estimated 

Error 

0.2 

0.6 

1 .17 

1 .10-10"4 

1 .16-1 0~4 

0.5 

0.2 

0.89 

1.13- 10-4 

1 .90- 10-4 

0.5 

0.6 

1 .08 

4.64- 1 O-4 

4.50- 10-4 

0.8 

0.2 

0.63 

1.24-10'3 

1.7M0"3 

0.8 

0.6 

0.88 

1 .48 -10" 3 

1 .38*10" 3 

0.9 

. - - - ■ 

0.3 

0.54 

4. M0’3 

8.9- 10-3 

Analogous  results  were  obtained  for  the  same  problem  with  3u/3x  = 0 at 
x = 0 and  unchanged  boundary  conditions  on  the  remaining  sides  of  the 
unit  square. 

4.4  Parabolic  partial  differential  equations 
Consider 

ut  = f(t,  x,  u,  ux,  uxx)  (t,  xKfl 

u(a,  t)  = f i ( t) , u(b,  t)  = f 2 ( t ) , t ? 0 ; u(x,  0)  = h(x),  a s x j b 

where  U - {a  < x < b,  t > 0} . Assume  that  f,  f ^ , f^  and  h are  such  that 

u € Cs(«). 

The  definition  of  spaces,  norms  and  operators  for  our  operator 
formalism  is  left  as  an  exercise  for  the  reader. 


f 
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We  will  use  the  method  of  lines  to  discretize  the  problem  in  space 
and  then  solve  the  system  of  ordinary  differential  equations  so  obtained 
with  two  simple  methods;  the  explicit  Euler  method  and  the  implicit  backward 
Euler  method.  It  is  plausible  that  for  sufficiently  small  h (discretization 
parameter)  there  exists  a smooth  expansion  of  the  global  discretization 
error  for  both  these  methods  (Stetter  (1965),  Keller  (1970)).  However, 
the  main  advantage  of  implicit  methods  is  that  one  can  use  large  time  steps 
and  still  obtain  accurate  results.  For  realistic  stepsizes  in  time  there 
are,  to  my  knowledge,  no  general  results  on  the  existence  of  smooth  expan- 
sions for  the  global  discretization  error.  The  system  of  differential 
equations  that  was  obtained  by  the  method  of  lines  is  stiff  and  hence 
results  on  error  expansions  for  numerical  methods  for  stiff  systems  of 
ordinary  differential  equations  may  be  useful  for  our  methods.  For  such 
problems  a discussion  of  error  expansions  for  large  values  of  the  stepsize 
h (not  only  asymptotically)  for  implicit  one-step  methods  can  be  found  in 
Dahlquist,  Lindberg  (1973).  These  questions  warrant  further  study  that  are 
outside  the  scope  of  this  report. 

Here  we  assume  the  existence  of  sufficiently  smooth  error  expansions 
for  the  methods  and  stepsizes  we  use  and  proceed  under  that  assumption  to 
estimate  the  global  discretization  error  with  our  algorithm.  Numerical 
results  indicate  that  the  assumption  is  reasonable. 

4.4.1  The  method  of  lines  with  Euler's  method 


Discretize  the  problem  according  to 


64 


xi  = a + i*h;  i=0,  1.....N;  h = b-  ^--a- 


tj  3 j’k;  j = 0,  1,  ; k = c«h  c $ 0.5 


uij  5 u(v  V 


U'J'1  " "'J  . f(t  x u ui*lj  ~ Vlj  Ui+1  j ' 2“ij  * . o 

k V Ui  j ’ 2h  1 


2h 


i 3 1(1)N  - 1,  j = 0,  1,  ... 


uiQ  - h(x.)  =0  i = 0(1 )N 

U0j  " = ° UN j W = ° J = 2*  ••• 

This  is  an  explicit  method  so  we  simply  proceed  one  time-level  a time 
computing  according  to  the  formula  above. 

The  perturbations  are  defined  according  to: 

Introduce  £ = U-j  i = 0(1)N,  j 3 0,  1,  ...) 

rl  >2  /,\  a j.  „ r nv^  >4  ( r\  ny2>4 


*uU)  = DT'*%  U)  - f(tr  X DX'’H  (0.  DX‘’H  (C)) 

ij  xi,tj  j i U Vlj  x1,tj 

i = 1 (1  )N  - 1,  j = 0,  1,  ... 


Here  the  operators  DT7’™,  DX7’™  are  of  the  type  discussed  in  section  3. 

X y L X j t 

For  the  interior  grid  points  we  have  e.g. 


DTx,2t  U) 
Xi’tj 


2k 


ox1.4  , , , ^1-2,1  - ! 8ti+l.i  ~ Ci*2j 

x.,t.v  ' 12h 

' *J 

The  perturbed  problem  is  solved  in  parallel  with  the  unperturbed  problem. 

Several  alternative  perturbations  exist,  cf.  section  4.1. 

The  problem  below  was  solved  on  an  IBM  360/75  using  double  precision 

2 

arithmetic.  We  used  N = 10,  and  c = 2.5/tt  . 
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ut  - uxx  0 < x < tt,  t > 0 u(0 , t)  = u(n,  t)  = 0 

with  the  following  initial  value  functions  h(x) 

(1)  h(x)  = sin(x) 

exact  solution:  u(x,  t)  = e't*sin(x) 

(2)  h(x)  = x(tt  - x) 

exact  solution:  u(x,  t)  * - Z (2n  - 1)~3  e"^n-l)  t Si n( 2n  - l)x 

" n=l 

(3)  fx  0 <;  x sn/2 

h(x)  = ( 

1 TT  " X Tl/ 2 < X < IT 

oo  . . 2 

exact  solution:  u(x,  t)  = - I (-l)^01-1'^2 — 1-  sin(m  x)-e”m  * 

" m=l(2)  r/ 

The  tables  below  give  the  maximal  errors  for  selected  time  levels  for  the 
different  problems. 


Problem 

#1 

N = 10,  k = 0.025 

t 

Actual  Error 

Estimated  Error 

0.05 

2.05-10-4 

2.13-10"4 

0.20 

7.05- 10-4 

7.46-10'4 

0.50 

1 .30-10'3 

1 .35-1 0"3 

0.75 

1 .52-10"3 

1 . 56* 10~3 

Problem  #2 

N = 10.  k = 0.025 

t 

Actual  Error 

Estimated  Error 

0.05 

2.6- 10-3 

8.4- 10-3 

0.20 

2.2- 10-3 

3.4- 10-3 

0.50 

3.2-10"3 

3.6- 10'3 

0.75 

3.9-10"3 

-3 

4.3-10 
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Problem 

#3 

N = 10,  k = 0.025 

~! 

t 

Actual  Error 

Estimated  Error 

0.05 

3.7-10'2 

1.2-10'2 

0.20 

1 .6-10'2 

4.9-10'3 

0.50 

9.8-10'3 

3. MO'3 

0.75 

8.3- 10-3 

3. MO'3 

For  these  problems  the  amount  of  work  needed  to  get  the  error  estimate  is 
approximately  equal  to  the  amount  of  work  needed  to  solve  the  unperturbed 
problem. 

From  the  numerical  results  it  is  obvious  that  the  estimates  are 
not  very  reliable  for  test  problem  3 where  the  initial  value  function  is 
not  smooth;  for  the  other  two  problems  the  error  estimates  are  quite  good 


4.4.2  The  method  of  lines  with  the  backward  Euler  method 
Discretize  the  problem  according  to 

xi  = a + ih;  i = 0(1 )N  h = (b  - a)/N 

tj  = j-k;  j = 0,  1 , .. . ; k = c-h 

uij  “ u(v  V 


' “ij  fit  X u ^ 

k 1 j+r  V uij+r 


UM  j+1  ' Vl  j+1 

2h 


Vl  j+1  ~ 2uij+l  * ui-l  j+l 


) = 0 


i = U1)N  - 1,  j = 1,  2,  ... 

ui0  ' h(xi)  = 0 1 = °<1)N 

u0j  - W * 0 uNj  - f2(tj>  = 0 


j = 1,  2,  ... 
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The  system  of  nonlinear  equations  that  is  obtained  in  each  time 
step  is  solved  by  Newton  iteration.  As  initial  guess  we  use  the  solution 
at  the  previous  time  level  and  the  iterations  are  terminated  when  the  last 
correction  is  less  than  e.  The  maximum  norm  is  used  to  measure  all 
quantities. 

Introduce  the  notation 

£ = (Cij;  i = 0(1 )N.  j = 0,  1,  ...) 

The  basic  discretization  for  the  perturbation  operator  is 


lj(t)  f(tj,  x.,  t...  DX'>«  <E),  <■«  (t)) 


1 J 


1 J 


i = 1 (1  )N  - 1 j = 1,  2,  ... 


v m 

The  operators  DX  ’ are  of  the  type  discussed  in  section  3. 

A } l 

For  the  perturbed  problem  we  use  modified  Newton  iteration  to  solve 
the  system  of  nonlinear  equations  obtained  in  each  time  step.  As  initial 
approximation  we  use  the  solution  of  the  unperturbed  problem  and  terminate 
the  iterations  when  an  estimate  of  the  iteration  error  is  either  less  than 
y times  the  estimated  discretization  error  at  the  current  time  level,  or 
less  than  e. 

The  problems  below  were  solved  in  this  fashion. 

(1)  ut  = (2  + sin(uxx))uxx  + (1  - sin(u) )u 

0 < x < "/2,  t > 0 

u(0,t)  = 0 ; u(V2,  t)  = e-t  t > 0 

u(x,  0)  = sin(x)  0 S x { »/2 
exact  solution:  u(x,  t)  = e_t  sin(x) 


(Kolar  (1972)) 
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(2,  3,  4)  Burger's  equation 

u = - u u + vu  0<x<tt,  t>0 

t x xx 

u(0,  t)  = u(n,  t)  = 0 

u(x,  0)  = sin(x) 


exact  solution: 


2 1 

i exp(-  v n t)n  In(^-) *sin(n  x) 
u( x , t)  = 4v — 1 

I0(2v')  + 2 E exP(‘  v p2  t)In(2v)*cos(n  x) 

n=l 

The  functions  I , n = 0,  1,  ...  are  the  modified  Bessel  functions.  (Cole 
(1951)).  We  have  used  the  following  values  of  v: 

(2)  v = 1 

(3)  v = 1/8 

(4)  v = 1/16 

In  the  calculations  for  the  tables  below  we  have  used  e = 10"6,  k = 0.1  and 
some  different  values  of  N.  The  computations  were  done  in  double  precision 
on  an  IBM  360/75.  The  maximum  errors  at  some  different  time  levels  and  the 
number  of  iterations  for  the  unperturbed/perturbed  problems  are  recorded  in 
the  tables. 


Problem  #1 , N = 10 , k = 0 . 1 


t 

Actual  Error 

Estimated  Error 

Number  of  iterations 

3.2-10"3 

2.6- 10"3 

5/2 

6.6.10"3 

5.8-10"3 

6/2 

1 .0 

4.0.10-3 

4.2-10"3 

6/2 

1.5 

2.22* 10-3 

2.19-1 0-3 

5/2 

2.0 

1 .22* 10"3 

1 . 21  - 10-3 

5/2 
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Problem  #4,  N = 20,  k =0-1 


t 

Actual  Error 

Estimated  Error 

Number  of  Iterations 

4.5*1 0- 3 

4. M0'3 

3/2 

0.5 

1 .8* 10'2 

1 . 7 • 1 O' 2 

3/2 

1 .0 

2.6-10'2 

2.3*1 0~2 

3/2 

2.0 

1.9-10'1 

3.5*1 0"^ 

3/2 

4.0 

1 .M0"1 

1.4-10'1 

3/2 

6.0 

3.9- 10-2 

3.8-10'2 

3/2 

8.0 

2. M0'2 

2.2-10'2 

3/2 

All  the  estimates  above  are  very  good,  except  for  problem  #4  where  for 
t ' 2.0  we  get  estimates  that  are  approximately  twice  the  actual  error. 

This  problem  developes  a fairly  sharp  front  with  large  derivatives  with 
respect  to  x.  When  the  front  disappears  the  estimates  become  quite  reliable 
again. 

Note  that  the  amount  of  work  needed  to  find  the  error  estimate  is 
quite  small  compared  to  the  amount  of  work  needed  to  find  tne  approximate 
solution. 

For  large  values  of  N we  need  approximately  N /3*  (number  of 
iterations  for  the  unperturbed  problem)  operations  to  find  the  approximate 
solution  and  N * (number  of  iterations  for  the  perturbed  problem)  operations 
to  get  the  error  estimate. 

4.5  Hyperbolic  partial  differential  equations 

For  the  commonly  used  discretizations  of  initial-boundary  value 
problems  for  hyperbolic  partial  differential  equations,  I know  of  no 
extensive  discussion  of  the  existence  of  smooth  expansions  of  the  global 
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discretization  error.  Some  very  interesting  results  are  given  in  Gourlay, 
Morris  (1968)  and  Skollermo  (1975a,  1975b).  In  essence  the  main  difficulty 
seems  to  be  how  to  represent  the  numerical  boundary  conditions  (i.e.  the 
extra  boundary  conditions  imposed  by  the  numerical  scheme)  so  the  errors 
due  to  that  representation  (these  errors  are  not  smooth)  is  sufficiently 
smal  1 . 


For  the  method  of  characteristics  some  authors  have  used  extrapolation 
to  increase  the  accuracy,  see  Lister  (1960),  Werner  (1968),  Smith  (1970). 

The  results  of  this  section  are  of  an  experimental  nature,  and  no 
rigorous  mathematical  analysis  is  attempted.  The  numerical  results  indicate 
that  by  careful  choice  of  the  perturbation  operators  one  can  obtain  good 
error  estimates  for  problems  with  smooth  solutions.  A further  study  of  this 
class  of  problems  may  prove  very  fruitful. 

Consider  the  initial-boundary  value  problem 
ut  = a ux  t 5 0,  0 s x s 1 
a > 0 

u(x,  0)  = f(x)  0 « x f 1 
u(l , t)  = h(t)  t i 0 
and  use  the  discretization 

x.  = i-h  i = 0,  1 , ... , N ; h = 1/N 

t j = j • k k = 0,  1 , . ..  ; k = c-h 


ui j+1  ' ui j-1  . Vl.i  " ui-l j . 

2k  ' a 2h 

uiQ  = f(xi)  i = 0(1 )N 


UNj  = h<tj> 


j - 0,  1,  ... 
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This  is  the  leap-frog  method.  Note  that  to  be  able  to  compute  the 
approximate  solution  with  the  formula  above  we  must,  in  some  way,  find 

i = 0(1 )N  (extra  starting  values) 


and 


uQj  j = 1,  2,  ...  (numerical  boundary  values) 

If  the  extra  starting  values  are  computed  with  a sufficiently  accurate 
method  we  have 


“u  = u(v  V * h‘ 


e(x.,  tj)  + hy 


1J 


+ R. 


1J 


where 


| R. . | = 0(hm)  m = min(4,  u + 1 ) 

* J 

u depends  on  the  method  we  use  to  find  the  numerical 
boundary  values 

e .j j = (-1)1  C(x^ , tj)  where  C is  a smooth  function 
(adapted  from  Skollermo  (1975a)). 

We  have  made  some  numerical  experiments  with  our  estimation  algorithm  for 
this  discretization. 

To  construct  the  perturbation  operator  <t>j^  we  proceed  in  the 
following  way: 

Defi ne 

5 = (f^.,  i = 0(1  )N,  j = 0,  1,  ...) 

e = (e-jj,  i = 0(  1 )N , j = 0,  1,  ...) 

6(U  = ((^i+lj  - Ci.lj)/2h,  i = 1 (1  )N  - 1,  j = 0,  1,  ...) 

Note  that  when  we  apply  6 to  the  non-smooth  error  term  hy*e  of  the  expansion 
above  we  get 
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6ij(hw  e)  = (-1)1’1 • C x ( X i , tj)-hw  + 0(h2+y) 

As  the  basic  discretization  for  the  perturbation  operator  ^ we  now  use 

« DTj’4  (?)  - a DX  (6(e)) 

i J x-  » x . > 


1 J 


1'  J 


where 


. 

DX  (6(e))  = E O.  6 .(e) 
xi,tj  s-4i  1S  SJ 


is  a linear  combination  of  some  of  the  components  of  6(e)  such  that 


DX  (6(a.  z))  ’s  consistent  of  order  4 with  z (x. , t.).  The  c 'fficients 

^ i * j ^ x i j 


a..  are  independent  of  h so 

hu  ^j(e)  = (-l)i'1(-ct(xi,  tj)  - a cx(x.,  t.))hy  + 0(hy+4) 

- 0(hu) 

This  choice  of  the  perturbation  operator  ^ insures  that  no  serious  loss 
of  accuracy  in  the  error  estimate  is  caused  by  the  irregular  error  term 
hy  • e . 


The  theorems  of  section  2 do  not  cover  this  kind  of  error  expansions 
and  perturbation  operators,  but  they  could  easily  be  modified  to  do  so. 

The  following  problem 

u . = u 0 < x < 1 , t > 0 

W»  A 

u(x,  0)  = sin(2w  x)  0 < x < 1 
u(l , t)  = s i n ( 2tt  t)  t > 0 

with  the  exact  solution  u(x,  t)  = sin(27i(x  + t))  was  solved  with  the  leap- 
frog method. 

First  we  extended  the  solution  to  the  half  plane  t z 0 and  used  the 
periodicity  of  the  exact  solution  to  find  the  numerical  boundary  values,  i.e. 
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U0j  = UNj 


1,  2,  ... 


We  also  used  the  periodicity  to  simplify  the  formulas  for  the  perturbation 
operator  <^. 

This  problem  was  included  to  see  if  our  algorithm  would  work  in 
the  case  where  the  numerical  boundary  condition  did  not  introduce  a non- 
smooth error  term.  Two  ways  of  finding  the  extra  starting  values  Up , 
i = 0(1  )N  were  tested, 

I.  Up  = sin(2n(x.  + k))  i = 0(1)N 

i.e.  values  from  the  exact  solution. 

H.  Up  = uiQ  + j x(ui+]  Q - u.p  0)  + \ X2(ui+1  0 - 2ui0  + ui_1  Q) 

with  X = k/h 

i.e.  the  Lax-Wendroff  scheme. 

In  tables  1 and  2 the  maximal  errors  for  some  time  levels  are  given  for 
these  cases. 

In  the  next  experiment  we  did  not  use  the  periodicity  of  the  exact 
solution,  but  used  the  following  two  formulas  to  find  the  numerical  boundary 
values 

A.  U0.M  ■ u0j  * i(utj  • u0j>  ; x ■ k/h 

i.e.  the  explicit  Euler  scheme. 

“•  u0j*l  * U1  j+1  • x(ulj*l  • W ' u0,j  + ulj 


- »(U,J  - u0J) 


X = k/h 


u . 


i.e.  the  Box  scheme. 

The  Lax-Wendroff  scheme  was  used  to  get  the  extra  starting  values 
p,  i = 0(1  )N.  In  tables  3 and  4 the  maximal  error  for  some  time  levels 


are  recorded. 
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In  all  the  calculations  for  the  tables  we  used 
h = 0.025  k = 0.01875 

and  the  computations  were  carried  out  in  double  precision  on  an  IBM  360/75, 

Table  1,  starting  procedure  I 


Estimated  Error 


1.4M0'3 
8.27- 10-3 
1.7M0'2 
3.42- 10-2 


Table  2,  starting  procedure  II 


Estimated  Error 


1.71-10  " 

2.18-10 

8.52- 10-3 

9.29-10 

O 

o 

1 

ro 

1 .71-10 

3*41- 10"2 

3.42-10 
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Table  4,  boundary  scheme  B 


t 

Actual  Error 

Estimated  Error 

0.15 

1 .7M0"3 

ro 

i 

o 

C\J 

0.75 

8.37-10'3 

9.09-10'3 

1.5 

1.69-10  e 

1 . 73-1 0-2 

3.0 

3.37-10"2 

3.45- 10-2 

The  nonlinear  problem 

ut  = fx  ^ u2)  0 < x < 1 , t > 0 

u(x,  0)=l-x  0<x<l 

u( 1 , t)  = 0 t > 0 

with  the  exact  solution 

u(x,  t)  = (1  - x)/(l  + t) 

(Gourlay,  Morris  (1968))  was  solved  with  the  leap-frog  method.  We  used  the 
Lax-Wendroff  starting  procedure  and  the  Box  scheme  for  calculation  of  the 
numerical  boundary  values. 

The  discretization  error  was  estimated  as  above.  For  h = 0.05 
and  k = 0.015  the  maximal  errors  for  some  time  levels  are  recorded  in  the 


table  below. 
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Table  5 


t 

Actual  Error 

Estimated  Error 

0.15 

3. MO'5 

3.3-1 O-5 

0.75 

3. MO'5 

3.6- 10"5 

1 .5 

2.2- 1 O-5 

3.2-1 O-5 

3.0 

1 .9- 10"5 

3.2- 1 O'5 

4.5 

2.0-10'5 

3.9-10'5 

6.0 

1 .9-1  O' 5 

4.6-1  O' 5 

7.5 

2.0*1  O'5 

6 . 0 - 1 O'5 

4.6  Integral  Equations 

First  consider  Fredholm's  integral  equations  of  the  second  kind 

b 

y(x)  - X / K( x , t)  y(t)  dt  - f(x)  = 0 
a 

with  the  following  discretization 

t1  * xi  = a + i.h,  i = 0(1 )N  ; h = (b  - a)/N 

N 

MO  = U.-'hx  I a.  K(x . , t.)  r.  - f(x.)  ; i = 0(1)N) 

M 1 j=0  J 1 J J « 

where  ctQ  = aN  = 1/2;  = 1 , j = 1 (1  )N  - 1 . 

We  have  used  the  trapezoidal  rule  to  approximate  the  integral  in  the  equation 
above. 

Construct  the  perturbation  operators 
N 

♦htVU)  = (Vhx  z0  ajv  K(xi’  tj)  Cj  - f(xi)  ; i =0(1  )N) 
v = 1,  2,  ... 
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Here  the  coefficients  a.  are  such  that 

J V 

N b s . +1 

z a.  ip(t.)  = / v(t)  dt  + z f.(*)  hJ  + 0(hs  ') 

j=0  JV  J a j=2( v+1 ) J 

The  system  of  linear  equations  obtained  for  the  unperturbed  problem 
is  solved  by  Gaussian  elimination.  The  LU  factorization  of  the  coefficient 
matrix  is  saved  for  use  in  the  '■•■'■'uence  of  perturbed  problems. 

The  unperturbed  problems  takes  approximately  (N  + 1)  /3  operations 
(ignoring  the  number  of  operations  needed  to  compute  K(x.,  t),  i = 0(1)N, 

* J 

j = 0(1  )N  and  f(xi),  i = 0(1)N)  while  each  of  the  iterations  in  the  iterative 

improvement  takes  approximately  2(N  + 1)  operations  (ignoring  the  number 

of  operations  needed  to  calculate  the  coefficients  a.  of  the  perturbation 

J * 

2 

operators,  approximately  3 [2(v  + 1)]  ). 

The  a-  are  best  computed  as  weights  of  a composite  quadrature 

J * 

formula  of  order  at  least  2(v  +1).  The  length  of  the  intervals  for  the 
basic  quadrature  formulas  are  M*h,  where  M has  to  be  a factor  of  N (M  may 
be  equal  to  N).  The  coefficients  of  the  basic  formulas  can  be  obtained 
as  the  solution  of  Van  der  Monde  systems  of  linear  equations  in  the  same 
way  we  obtained  the  coefficients  of  the  differentiation  formulas  of  section 
3. 

In  Van  der  Sluis  (1972)  a discussion  of  asymptotic  error  expansions 
for  quadrature  formulas  of  this  kind  can  be  found,  and  in  Laurent  (1964), 
Stetter  (1965)  asymptotic  expansions  for  the  global  discretization  error 
for  our  method  are  discussed. 

The  problems  below  were  solved  on  an  IBM  360/75  in  double  precision. 
All  problems  are  taken  from  Netravali,  de  Figueiredo  (1974). 
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f 

1 

(1)  y(x)  - / x t y(t)  dt  - (ex  - 1)  = 0 

0 


exact  solution:  y(x)  = e 

1 

(2)  y(x)  - / x t y(t)  dt  - (sin(nx)  - x/tt)  = 0 

0 


exact  solution:  y(x)  = sin(nx) 

(3)  y(x)  - / x4  ext  y(t)  dt  - (x  - x3[ex(l  - 7)  + 7])  = 0 

0 xx 


exact  solution:  y(x)  = x 

1 Ax 

(4)  y(x)  - / x4  eXt  y(t)  dt  ~ (sin(irx)  - ie  +y  h = 0 

0 + IT 


exact  solution:  y(x)  = sin(iTx) 

We  have  used  N = 16  and  record  below  the  maximum  errors  in  the  successive 


i terates . 


Problem 

Error  in  Iterate  Number 

0 

1 

2 

3 

4 

5 

6 

1 

2.1 

IQ'6 

-9 

2.1  10  * 

-12 

2.0  lU 

1.6  10'15 

-15 

7.1  10  13 

2 

1.4 

IQ'6 

2.7  10'11 

2.7  10‘14 

4.8  10‘15 

3 

2.3 

IQ’5 

1 .8  10"9 

-13 

1.6  10  1 

-15 

7.5  10 

4 

6.6  10'3 

5.7 

IQ'5 

4.8  10"9 

4.3  10"11 

-13 

3.7  10  1-3 

Note  1 The  coefficients  a,..,  j = 0,  1,  N,  v = 1,  2,  v can 

JV  maX 

also  be  made  independent  of  v if  they  are  chosen  as  a.  = w.,  j = 0, 

J v j 


where 
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N b 

Z w *(t  ) = f i,(t)  dt  + 0(hy) 
j=0  J J a 

and  u >.  (vmax  + l)-2. 

Note  2 The  technique  described  above  can  be  used  for  other  discretizations 
of  the  integral  equation,  using  e.g.  non-uniform  grid  points  , i =0,  1, 
....  N.  We  can  also  use  the  same  technique  for  other  types  of  integral 
equations,  e.g.  nonlinear,  as  long  as  the  assumptions  of  theorem  3 or  4 are 
satisfied.  Essentially  we  require  smoothness  of  the  exact  solution  and  a 
smooth  error  expansion  for  the  solution  of  the  discretized  problem. 

Consider  nonlinear  integral  equations  of  the  type 
b 

y(x)  - f K(x,  t,  y(x) , y( t) ) dt  - f(x)  = 0 
a 

with  the  following  discretization 

t.  = x.  = a + i*h  , i = 0(1)N  , h = (b  - a)/N 

N 

*hU)  = - h dj  K(Xi  , ty  Cj)  - f(x.)  ; i = 0(1)N) 

where 

a0  = aN  = 1/2  ; OL.  = 1 , j = 1(1  )N  - 1 

The  system  of  nonlinear  equations  obtained  is  solved  by  Newton  iteration. 

The  perturbations  are  constructed  analogously  to  the  perturbations 
for  the  linear  problems  above,  and  the  system  of  nonlinear  equations  obtained 
for  each  perturbed  problem  is  solved  by  Newton  iteration.  The  iterations 
are  terminated  when  the  iteration  error  becomes  less  than  e. 

The  problems  below  were  solved  in  this  fashion  on  an  IBM  360/75  in 


double  precision. 
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(1)  y(x)  - / x-t-y(t)2  dt  - 3x/4  = 0 

0 

exact  solutions:  y(x)  = x and  y(x)  = 3x 

(Moore  (1968)) 

1 2 

(2)  y(x)  - f y(t)2  dt  - ex  - e2x(e2  - 1 )/ (4(1  + x))  = 0 

0 

exact  solution:  y(x)  = ex 

(arti ficial ) 

-14 

We  have  used  N = 10,  e = 10  and  recorded  below  the  maximum  errors 
in  the  successive  iterates. 


Problem 

error 

in  iterate  number 

0 

1 

2 

3 

4 

5 

6 

1 

0.5-10'6 

0.5-10  0 

0.5-lCf 10 

-12 

0.5-10 

0.5-10"14 

2 

3.0* 10-6 

4.3-10'8 

9.3-10'11 

8.7-10'11 

Note  1 With  N = 10  one  cannot  expect  to  get  the  error  less  than  OJh11), 
which  is  obtained  after  five  iterative  improvements.  For  problem  1,  however 
the  sixth  iterate  improves  the  accuracy  as  much  as  the  previous  iterates. 
This  astonishing  result  is  due  to  the  fact  that  in  the  error  expansion 

v(a.  y)  = A°{F(y)  + i hj  *.(y)}  + 0(hM+1) 
n,V  n n .Hv+1)-2  j 

we  have  ^.(y)  = 0 for  the  exact  solution  y of  the  integral  equation.  This 

J 

happens  because  the  quadrature  rules  that  are  used  for  v = 1,  2,  ...  are 

exact  for  polynomials  of  degree  three  or  less  and  for  the  exact  solution 
3 

the  integrand  is  t . cf.  note  1 after  the  numerical  results  of  section 
4.3.1. 
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5 . Concluding  remarks 

In  the  final  stage  of  the  work  with  this  report  we  got  a preliminary 
version  of  a paper,  "Iterated  defect  corrections  based  on  estimates  of  the 
local  discretization  error,"  by  R.  Frank,  J.  Hertling  and  C.  W.  Uberhuber, 
report  number  18/76  from  the  Institut  fur  Numerische  Mathematik,  Technical 
University  of  Wienna,  Wienna,  Austria.  There  the  authors  consider  an 
algorithm  very  similar  to  our  algorithm  for  iterative  improvement.  However 
their  perturbation  operators  (or  approximations  to  the  local  discretization 
error)  are  computed  differently.  For  the  two  point  boundary  value  problem 
y"  = f(t,  y)  y (0)  = A y(l ) = B 
with  the  discretization 


Eo  - A 


♦h(E)  ■ 


Ef*I 


2Ei  * Ei-1 


- fit,.  V 


(•1,2 N - 1 


\E»  - B / 
they  define  locally  smooth  functions  P^(t,  n^*) , k = 1,  2,  ...,  N - 1 that 

interpolate  the  solution  of  <t>.  ( n® ) = 0 at  some  points  surrounding 


x = . Then  they  compute  the  local  discretization  error  as 

V"0)  ■ <pX-r  "°>  - 2P°<V  p°>  p pk‘Vc  "0»'h2  - <P°>"<V  p°> 


k = 1,  2,  ....  N - 1 


The  succeeding  approximations  to  the  local  discretization  error  (corresponding 
to  our  approximations  4^  v(nv  '),  v = 2,  3,  ...)are  computed  by  the  same 
technique,  but  now  higher  and  higher  order  interpolation  is  used  to  define 
the  functions  P^( t , n).  The  polynomials  P^( t , nV)  do  not  need  to  be 
computed  explicitly  but  the  second  derivative  can  be  computed  by  forming 
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linear  combinations  of  the  components  of  n , as  we  do  when  we  approximate 
linear  functionals. 

In  our  notation  (cf.  section  4.2)  they  define  families  of  perturbation 
operators  <f>,  , v = 1 , 2,  . . . according  to 

i"'\  \ 

Vvtt)  'I  tk~1  V Vl  - °t,2t2V(^)  k - 1.  2 N - 1 


'N-l  ' 6 


Note  that  for  any  z € E 


(a.  z)  = a“{F(z)  + i f.(z)  h2j}  + 0(hM+1) 


where 


V2>  - 


,(2j+2) 


i + 2 ! 


♦f,  v(Ah  z)  = Z f.(z)  h2j  + E f - ( z ) h2j}  + 0(hM+1) 
,V  n j-1  3 j=v+l  VJ 


Vj-1  'j’ 


(if  D2’2  2v  uses  points  symmetrically  distributed  around  t,)  with  proper 
definitions  of  A.  and  a?.  The  improved  solutions  are  computed  according 


4>h(  n ) = 0 


*h(r,1)  ‘ Vi^"^  = 0 i=1.2. 
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Our  theory  does  not  cover  this  type  of  perturbations,  but  theorems  similar 
to  ours  could  be  proved  with  our  technique  for  this  algorithm. 
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